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FOREWORD 


This report contains a description of a computer code for solving two 
hypersonic vehicle trajectory optimization problems, that is, maximum cross¬ 
range and maximum plane change. The work was performed by David G. Hull 
and Jason L. Speyer of the Department of Aerospace Engineering and 
Engineering Mechanics, University of Texas at Austin, Austin, Texas, 78712. 
The program was sponsored by the Air Force Flight Dynamics Laboratory 
under Contract F33615-79-C-3030 and work unit number 2404 07 32 
and was managed by Dr. L. Earl Miller. 
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SECTION I 


INTRODUCTION 

The purpose of this study is to certify the ability of a particular 
timizatiop code to solve trajectory optimization problems associated with 
hypervelocity vehicles, that is, the maximum crossrange problem and the 
maximum plane change problem. The optimization problems are converted into 
parameter optimization problems (nonlinear programming problems) and are 
solved by the augmented-Lagranglan method. The particular optimization code 
which will be used has been created at the Atomic Energy Research 
Establishment in Harwell, England. 

The physical model used for the trajectory problems is described in 
Section II. In Section III, the optimal control problem is converted into 
a parameter optimization problem, and the optimization code is discussed 
in Section IV. That part of the code to be provided by the user is 
presented in Sec'Ion V. The solution of the optimization problem is 
carried out in Section VI, and Che conclusions are presented in 
Section VII. Finally, Che documentation provided with the optimization 
code, Che listing of the optimization code and user code, and Che 
definition of the variables in the user code is presented in the 
appendices. 
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SECTION II 


PHYSICAL HODEL 

In this section, the physical model used for the reentry and plane change 
problems is discussed. The model includes the equations of motion, the earth, 
the atmosphere, the aerodynamics, the propulsion, the performance Indices, the 
boundary conditions, and the inequality constraints. 

2.1 Equations of Motion 

The equations of motion used for the reentry and plane change problems are 
those for thrusting flight over a rotating earth. If the vehicle is assumed to 
be flying west to east and if a positive bank angle is assumed to generate a 
heading toward the north, these equations are given by 

8 ■ V cos y COB i///r cos f 

(J ■ V cos Y sin 'I'/r 

r ■ V sin Y (1) 

V ■ (T cos a - D)/ra - g sin y + w^r coo ^ (sin y cos ^ - cos y sin ^ sin i/;) 

Y “ (T sin a + L) cos p/mV + (V^/r - g) cos y/V + 2u cos ^ cos ij; 

+ (u^r/V) cos (cos Y cos ^ + sin y sin ^ cos ti>) 

■ (T sin a + L) sin p/mV cos y - (V/r) cos y cob ij/ tan ^ 

+ 2a) (tan y cos ^ sin ifi - sin - (a)^r/V cos y) sin ^ cos | cos 4/ 

In these equations, e is the longitude, (P is the latitude, r is the radial 
distance from the center of the earth to the vehicle center of gravity, V is 
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the velocity relative to the earth, y is the flight path Inclination, 
is the heading angle, m is the mass, T is the thrust, D is the drag, L Is 
the lift, a l8 the angle of attack, g is the bank angle, and u is the 
angular velocity of the earth. 

2.2 Earth 

The earth Is assumed to be a sphere whose radius represents mean sea level 
and Is denoted by rg . As a consequence, the acceleration of gravity is given 
by the inverse-square law 


g - gj.(rg/r )2 ( 2 ) 

where gg is the acceleration of gravity at sea level. Also the altitude of 
the vehicle above mean sea level satisfies the relation 

h • r - tg (3) 

The values of the constants associated with the characteristics of the earth 
are listed below: 

Tg - 20926428 ft 

gg - 32.174 ft/scc2 (4) 

u ■ 7.2921151E-5 rad/sec 

2.3 Atmosphere 

In order to obtain atmospheric properties as a function of altitude, the 
approximate atmosphere known as Che 1962 U.S. Standard Atmosphere has been 
employed, with the additional assumption that the composition of the atmosphere 
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is constant. Here, the atmosphere Is divided Into a nusiber of layers In 
which the temperature gradient is assumed constant. Hence, the absolute 
temperature In each layer Is given by 

T - T + 1 (h - h^) . h ^ h^ (5) 

n n n — n 

where Is the absolute temperature at the beginning of the layer, 1^ 

Is the constant temperature gradient, and h is the altitude at the beginning 

n 

of Che layer. Once the temperature Is know, the speed of sound Is obtained 
from the relation 

a - (kR (6) 


where k Is the ratio of specific heats of air and R Is the gas constant of 
air at sea level. For those layers where the temperature gradient Is nonzero, 
the density ratio, a - p/p^ , can be expressed as 


0 


C 

n 


-(1 + gc/R 1 ) 
T bn 


(7) 


where Pc Is the density at sea level and C Is a known constant for each 
^ n 

layer. For chose layers where the temperature gradient is zero, the density 
ratio becomes 


o • C exp (-g„h/R T ) (8) 

n on 

The values of Che constants associated with the atmosphere are presented 
below and In Table 1: 

k - 1.4 

R - 3086.9629 ftVsec* ’K 

pg - 2.3769E-3 slugs/ft* (9) 
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Table. 1. Atmospheric Constants 


Layer 

(n) 

(ft) 

"n 

(*K) 

^n 

(*K/ft) 

C 

n 

1 

0 

288.15 

-1.9812E-3 

3.401824655257E-11 

2 

36,089 

216.65 

0 

1.683376997149E+00 

3 

65,617 

216.65 

3.0480E-4 

9.8178589l4969E-f80 

4 

104,987 

228.65 

8.5344E-4 

1.506414967722E+29 

5 

154,199 

270.65 

0 

4.394749884481E-01 

6 

170,604 

270.65 

-6.0960E-4 

4.717851690435E-43 

n 

/ 

200,131 

252.65 

-1.2192E-3 

1.562793740651E-22 

8 

259,186 

180.65 

0 

5.028946984109E+01 


During the optlmlz.itlon process, trajectories can be obtained which go to 
very high or very low altitudes. To avoid a fatal error associated with expo¬ 
nential overflow, the following additional features have been Incorporated 
In the model atmosphere: 


h 1 0 : T - 288.15 - 1.9812E-3 h 

o - 1. - 2.9262913E-5 h 

h ^ 750,000: t - 180.65 

0 - 8.111816E-18 


( 10 ) 


Hence, below sea level, the density ratio Is assumed to vary linearly with 
altitude, and at high altitudes, it is assumed constant at the value for 


750,000 ft. 
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2.4 Aerodynamics 


The drag and the lift are related to the drag coefficient and the 
lift coefficient C as follows: 


D - (1/2) OSV^C, 


L - (1/2) pSV^Cj^ 

where S Is the aerodynamic reference area. The drag and lift coefficients 
can be expressed in terms of the axial and normal force of coefficients as 


C_ ■ C. cos a + C„ sin a 
DA N ■ 


COB a - sin o 

The axial force coefficient is composed of a skln-frictlon term and a pressure 
term, that Is, 


c* ■ c. ■*- C. 

^ ^SF \r 


where 


C, - A.h^ + B,h + C, 

Asf 1 11 


“ ^ 2 “ » 2 “ ^2 


In these relations, A^ , , and are known functions of the Mach 

number, M ■ V/a , while A 2 , and are known functions of Mach 

number and angle of attack. Finally, the normal force coefficient is given 
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S “ ^3“^ ®3“ ''' ^3 


(15) 


where » 83. and are known functions of Mach number. 

The values of the constants associated with the aerodynamics are given 
below and in Tables 2 throug!i 4. 


S - 125.84 ft? (16) 


Table 2. Axial Skin-Friction Force Constants 


M 


«1 

"1 

0.2 

0 

7.80F.-3 

.0114 

1.2 

0 

7.70E-8 

.0076 

5.0 

0 

5.60E-8 

.0028 

10. 

2.40E-12 

-7.11E-7 

.0578 

>. 20. 

1.38E-12 

-3.81E-7 

.0297 



Table 3. Axial 

Pressure Force Constants 



M 

^2 


®2 


^2 


a < 16* 

0 > 16” 

a < 16” 

a > 16' 

a < 16' 

a" > 16* 

0.2 

-l.OOE-4 

-l.OOE-4 

-2.19E-3 

-2.19E-3 

.0350 

.0350 

1.2 

-7.81E-5 

-7.81E-5 

-1.25E-3 

-1.25E-3 

.0760 

.0760 

5.0 

1.09E-5 

4.86E-6 

-9.62E-4 

-1.13E-4 

.0324 

.0204 

10. 

1.41E-5 

-6.60E-6 

-9.00E-4 

3.49E-4 

.0261 

.0114 

> 20. 

1.33E-5 

-6.25E-6 

-8.19E-4 

3.58E-4 

.0243 

.0105 









7 



Table 4. Normal Force Constants 


M 


»3 

s 

0.2 

1.33E-4 

2.68E-2 

-.030 

1.2 

-7.81E»-5 

2.90E-2 

-.030 

5.0 

2.94E-4 

1.41E-2 

-.035 

10. 

4.38E-4 

6.50E-2 

-.035 

^ 20. 

4.38E-4 

6.50E-2 

-.035 


In the computer program, the above tables are read by linear Interpolation. 
Also, for M 20 , the M •• 20 value Is used because the extrapolation of the 
above numbers could lead to wrong results. 

2.5 Propulsion 

For the plane change problem, rocket engines are used to aid In the plane 
change and to Increase the energy of the vehicle. It Is assumed that the engines 
bum once and that the propellant mass flow rate and the thrust are constant 
during the bum. Hence, during the bum, the mass of the vehicle satisfies 
the relation 


m - - 8(t - tj^) 


(17) 


where 6 Is the mass flow rate and where and t^ are the mass and the 

time at the initiation of the bum. 

The values of the constants associated with the propulsion characteristics 
are as follows: 
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nij^ ■ oIq •• 335.67477 slugs 
6 ■ 0.34768573 slugs/sec 

T - 3300 lb (18) 


n ■ 179.18195 slugs 
P 

At * 515.35606 sec 
o 

where m Is the propellant mass and At. Is the total bum time. In the 
p b 

reentry program, the propellant mass is assumed to have been expended, so the 
Initial mass Is m^ * 156.49282 slugs. 

2.6 Performance Indices 

In the reentry problem, it is desired to maximize the crossrange. Hence, 
to have a minimization problem, the performance Index is written as 


J - - (19) 

For the orbital plane change problem. It is desired to maximize the 
plane change. Maximizing the plane change Is equivalent to maximizing the 
orbit Inclination because the Initial orbit Inclination Is zero. The orbit 
Inclination Is the angle between the Inertial angular momentum vector 


H 


T X V 

^ inertial 


( 20 ) 


and the angular velocity vector of the earth, where the bar denotes a 
vector. From the definition of the dot product. It Is seen that 


cos 1 ■ (H ' (i))/Ha) 


( 21 ) 
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AIho, the Inertial velocity Is given by 


'^inertial 


V + « X r 


( 22 ) 


Finally, If the vector operations are carried out, the problem of maximizing 1 
becomes that of minimizing the cosine of the Inclination at the final point: 


J « (cos 1)^ 

where 


(23) 


cos 1 


cos ^ (V cos Y cos rv> cos ♦) 


[(V cos Y cos + riD cos ♦)^ + (V cos V sin 4')^] 


172 


(24) 


2.7 Prescribed Boundary Conditions 

The Initial conditions for both problems are as follows: 

tQ - 0 , % - 0 , “ 0 , h^ • 364,800 

Vq - 24,500 ft/sec , Yq - -1.2 deg , 'I'q “ 0 

For the reentry problem, the final condition is given by 

h^ - 100,000 ft 

The final conditions for the plane change problem are the following: 
h^ ■ 608,000 ft , 23,500 ft/sec , Yj ■ 0 


(25) 


(26) 


(27) 
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2,8 Inequality Constraints 


The angle of attack and the load factor, n * L/W , are required to 
satisfy the Inequality constraints 

a ^ 40 deg , n ■ £ 4,5 (28) 


2.9 Nondlmenelonal Variables 

Nondlmenslonallzatlon has the effect of scaling the variables and can be 
very beneficial to optimization. In this system, all angles are In radians. 
Furthermore, to simplify the notation, a nondlmenelonal variable is denoted 
by the same symbol as the dimensional variable and a tilde. Hence, the 
following nondimensIona1 variables are defined: 

(29) 

= g 


t/(rg/gg)^^^ ■ ^ ^ ^ ^ ^ “ 

T/mogg = T . D/n^gg -= D , L/m^gg = L , m/m^ h m , e(rg/gg) /m^ 


M 
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SECTION III 


FORMULATION OF THR OPTIMIZATION PROBLEM 

In chis section, the optimal control problem is stated In terms of the 
nondlmensional variables, and the parameter optimization problem Is formulated. 

3.1 Optimal Control Problem 

The optimal control problem Is to find the angle of attack and roll 
angle histories which maximize the crossrange for the reentry problem 
which maximize the orbit inclination (or minimize its cosine) for the plane 
change problem. In terms of the nondlmensional variables presented in 
Equation (29), the optimization problem can be stated as the following 
minimization problem: 

Find the control variable histories a{c) and v)(t) which minimize 
Che performance index 



subject to the differential constraints 

de/dt ■ V cos y cos ^/{r cos $) 
d$/dt ■ V cos y sin ip/r 
dr/dl • V sin y 

dV/dl - (T cos a - 6)/i-sln y/i^ + w^r cos $(Bln y cos $ 
- cos y sin ^ sin ip) 
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dV/dl ■ (T sin d + L) cos u/mV + V cos y/r ~ coi Y/r^V + 2w cos $ cos i 
+ (lu^r/V) cos $(co8 Y cos $ + sin y sin ♦ sin j)) 

(31) 

dij/dt ■ (T sin a + L) sin ;j/inV cos y “ (V cos y/r) cos il> tan ^ 

-t- 2ai(tan y cos ^ sin tji - sin 
- (Z^r/^ cos y) Bin ^ cos ^ cos ^ , 


subject to the prescribed boundary conditions 


t. = 0 . G. - 0 , “ 0 , r- - 1.017432502 


V - 0.94420438 . y^ - -2.0943951E-2 , k - 0 


Reentry: 


= 1.0047786464 


Plane Change : = 1.029054170 , * 0.90566542 , Y^ = 0 


and subject to the inequality constraints 


a _< 0.69813170 , n £ 4.5. 


For the reentry problem, the vehicle is gliding with the engine off. Hence, 
the thrust, the mass flow rate, and the mass satisfy the relations 


T«0, 6"0, m"l 


while the drag and lift are given by 


D - 15268.635 C..aV^ 


t - 15268.635 C oV^ 
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For the plane change prohlei&t the engines are Ignited at , whose optimal 

value is to be determined, and burn for the time period At, ■ 0.63901697 . 

b 

Here, the thrust, the mass flow rate, and the mass are given by 



e - 


0 engine off 

10.83533982 engine on 


m 


1. - S(t - f^) , 

0.46620367 


0 < t < t,. 


Ko 


(36) 


while the drag and the lift become 


D - 932.34367 
I • 932.34367 


(37) 


The reason for the different constants In Equations (36) and (37) is that the 
initial mass for the reentry problem is less than the Initial mass for the 
plane change problem (see Section 2.5). 

3.2 Parameter Optimization Problem 

In general, the optimal control problem of the previous section can be 
stated In standard matrix notation (Reference 1} as follows: 

Mlmlmize the performance index 

J - <t(tj, Xj) (38) 
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gubject Co Che dlfferentlsL constraints 


i, - f(t, X, u) , 

Che prescribed boundary conditions 

Cq - 0 , Xq i given , S»(tj, x^) - 0 , 
the control variable Inequality constraints 


( 39 ) 


(40) 


C(t, X, u) ^ 0 , (41) 

and the state variable inequality constraints 

S(t, x) > 0 . (42) 

In Che parameter optimization method which will be used to solve these 
problems, the inequality co . 'alnts must be in an algebraic form. This can 
be accomplished by forming an integral over the region where Che inequality 
constraint is vlolsted. In other words, if e(t) ^0 is a scalar inequality 
constraint, an appropriate Integral constraint is 

r^f 9 

E ■ - / (min (e, 0)]^ dt ^0 (43> 

This integral constraint can be converted into a differencial equation and an 
algebraic inequality constraint. Let 

\+l ' * ^n+1,0 " ° 

so chat Equation (43) becomes 


E 


X ., £ > 0 

n+l,f — 


(45) 


15 



If the Inequality constraints are converted into differential equations 
and algebraic constraintSt the optimal control problems can be restated as 
follows: 


Minimize the performance index 

J ■ ♦ (t^, y^) 

subject to 

y - f(t, y, u) 

Cq “ 0 . Vq = given, 


(46) 


(47) 


(tf. yp - 0 , (48) 

and to 

<3(t^, Yf) i 0 . (49) 

The vector y now contains the original state x and the states Introduced 
by converting the inequality constraints. 

The conversion of this optimal control problem into a parameter optimization 
or nonlinear progranmlng problem is performed in two steps. The first step is 
to convert the final time into a parameter, and the second step Ic to para¬ 
meterize the control histories. 

The final time can be made into a parameter by introducing the transformation 

T - t/t^ . (50) 
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This Cransfonnatlon allows the optimal control problem to be rewritten as 
Minimize the performance index 

( 51 ) 

subject to 


to 


y' - g(T, y, u, tj) 


(52) 


4' (y£, tj) - 0 , (53) 

and to 

0(y£. tj) i 0 . (54) 

Note that the prime denotes a derivative with respect to t and that this 
transformation fixes the Interval of integration. 

The control histories can be parameterized In a number of ways, but the 
simplest is to use a set of nodal points and create the function by linear 
interpolation (see Figure 1 where a and u are control variables). This 
approach allows the control to be written in the form 

Uj^ - u^(t, Sj^) (55) 

where a^ is a vector of parameters, that is, the nodal points defining the 
control history . 

At this point, the unknown parameters are lumped into a single vector 
X defined as 
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X ■ ^1* • * • • ®pl (56) 

where p is the numbir of control functions. For the plane change problem, 

X also contains the Ignition time, t^ , as the last parameter. The total 
number of parameters is N . If values are given to the components of X , 
that is, if the final time, the control histories, and the ignition time are 
knovm, the differential equations (52) cen be integrated from ~ ^ 

Tj - 1 to form the function ■ y£(X) . In turn, can be eliminated 

from the performance index and the constraints to form the parameter version 
of Che optimal control problem: 

Minimize the performance index 

J - P(X) (57) 

subject to the equality constraints 

C^(X) - 0 , i - 1 K (58) 

and the inequality constraints 

C^(X) >0, i-K + l-^M (59) 
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SECTION IV 


PARA>fETER OPTIMIZATION METHOD 


In this section, the theoretical heals of the method which has been 
selected for solving the optlrlzatlon problem Is discussed. Also, some 
remarks are made about the code. 


4.1 Augmented-Lagranglan Method 

The augmented>Lagranglan method is a particular formulation of the penalty 
function method which allows convergence to be achieved without having to drive 
the weights to Infinity (see Reference 2). Here, the performance Index Is 
defined as 

m 

J - PrX) + (1/2) I 0. l2(X, 0.) (60) 

1-1 ^ ^ ^ 


where 


c^(X) - 


, 1 “ 1 K 


r^(x, e^) 


(61) 


min [C^(X) - 0^, 01 , 1 K + 1 -► M 


and where Is a positive, constant weight. 

A description of the computational algorithm la as follows: 


1. Guess Initial values for X , ' Since x represents 

the final time, the controls, and perhaps tha Ignition time, it should 
be easy to run soma constant control trajectories and find a value for 
X which gives a low value to the performance index P . For the 
first run, the constants are usually oet equal to zero. Finally, 

the weights o^ can be obtained by requiring each constraint term In 
(60) to have the same order of magnitude as F . 
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2. Minimize the performance Index with respect to X holding and 

o^ constant. This step represents an unconstrained minimization 
which is performed with the Davldon-Fletcher-Powell method. 

3. Assuming that X - X(9) , minimize the perfomance Index with res- 

I 

pect to 0 , starting from the X which results from step (2). 

The purpose of this step Is to find new values of 6 such that the 
constraints are satisfied. Only a single iteration of the optimization 
Is performed. For equality constraints, this Is done with the 
Newton-Raphson method, while for inequality constraints, a minimi¬ 
zation problem Is formed which guarantees that the corresponding 
Lagrange multipliers remain nonnegative. 

4. If the rate at which a constraint is being reduced is not sufficiently 
large, the value of the corresponding weight 3^ Is increased. 

5. If convergence has not been achieved, go to 2. 

4,2 Optimization Con e 

The code which h»s been selected for solving the optimization problem has 
been obta. ned from the .'VVRE in Harwell, England. It is a general purpose code 
for solving the nonlinear programming problem using the augmented-Lagranglan 
method. The code contains five subroutines which are described briefly below. 

VFOIA: This subroutine directs the optimization process. 

VA09A: This subroutine performs unconstrained minimization using the 

Oavldon-Fletcher-Powell method. The metric is represented In a 
a factored form to preserve positive definiteness so that a crude 
one-dlmenslonal search can be used. Convergence is superlinear. 
NC13.A: This set of subroutines (HCllA, B, C, D, and E) performs matrix 
manipulations such as factorization and the metric update. 
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\'E04A: This subroutine performs the optimization for satisfying the 
constraints. 

VFOlZ: This subroutine fcrms the augmented Lagrangian (60) and Its 

derivative from the performance Index (57), the eqitallty constraints 
(58), the inequality constraints (59), and their derivatives. 

The documentation which has been supplied with these subroutines Is contained 
In Appendix A. Also, a listing of the code is presented In Appendix B. That 
part of the code which must be supplied by the user is shown in Appendix C. 

The user must provide a subroutine entitled VFOIB which computes the 
performance Index F and the constraints C as well as the derivatives Fy 
and . In order to minimize the amount of set-up time, it has been 
decided to compute the derivatives numerically, using central differences. 

Thia approach also allows the user to change the physical modal without having 
Co change derivative subroutines as well. 

This optimization code has been verified by applying it to the solution 
of a number of problems whose solutions are known. These problems include a 
simple quadratic function in two variables with a linear equality or inequality 
constraint and the Rosenbrock function with a quadratic constraint. Also, the 
optimal control problem known as the lunar launch problem has been solved. 

The problem is to transfer a rocket from the surface of Che moon to orbital 
conditions in minimum time, assuming that the ratio of the thrust to Che mass 
la constant. The control variable la the angle between the thrust vector and 
the horizontal. To check out inequality constraints, an active upper limit 
has been imposed on the steering angle. In all cases, the Harwell code per¬ 
formed extremely well, and the known extremal solutions were obtained in a 
correct or reasonable number of ICerationa. 
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After solving the above problems, it became apparent that some minor modi¬ 
fications had to be made in the Harwell code. First, if an inequality constraint 
w.ns Imposed and was exactly equal to zero on the first iteration, the progfAm 
aborted. To prevent this, the constraint was set equal to a small positive 
value. Second, for solving more complex problems, it became apparent that it 
would be necessary to interrupt the computation, store soma results on a tape, 
and restart the computation. Such modifications have been made. Finally, for 
solving the reentry and plane change problems, some procedure for preventing 
the flight path inclination from becoming -90 deg had to be Incorporated. 

If this happens, the heading angle equation blows up, and the program 
terminates. This situation occurs during the one-dimensional search associated 
with the unconstrained minimization. If too large a change is made in the 
parameters, the vehicle can end up in a vertical dive. The problem has been 
eliminated by monitoring the flight path inclination on every integration step. 

If it gets lower than -60 deg , the one-dimensional search is terminated, the 
search stepsize is reduced, and the search is restarted. 
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SECTION V 


USER CODE 

The user of the optimization code must provide a main program and a 
subroutine or sequence of subroutines In which the performance Index, the 
constraints, and their derivatives are calculated. These parts of the computer 
program are discussed here. Because of the similarity of the two trajectory 
optimization problems, It has been possible to combine them Into a single 
program. The listing of the user code is presented In Appendix C. The 
variables In the user code are defined In Appendix D. 

5.1 Main Program 

The main program MRRV performs two functions. First, all of the 
parameters needed by VFOIA are defined (see Appendix A). Second, the 
process for Interrupting the computation and storing Intermediate numbers 
on tape as well as reading Intermediate numbers from tape and restarting the 
computation Is controlled. 

5.2 Performance Index. Constraints, and Derivatives 

Subroutine VFOIB calls subroutine SG, which computes the performance 
index and the constraints, and subroutine SGX, which computes the derivatives 
of the performance index and the constraints. 

In subroutine SG, the tables defining the control variable histories 
are formed from the current values of the parameters. The Initial conditions 
for the differential equations are defined, and the differential equations of 
motion are Integrated to the final time. From the final values of Che states, 
the performance Index and the constraint residuals are computed. While the 
trajectory Is computed In nondlmenslonal variables, the printout Is In terms 
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of dimensional variables. For Che plane change problem, the engine is ignited 

at t. and shut off at t. + At. . Finally, if the flight path Inclination 
b b D 

becomes less Chan -60 deg at any point of the trajectory, the integration 
is stopped, and the computation is redirected to the one-dimensional search 
for a decrease in the search stepsize. 

The integration is performed using a fixed-step, fourth-order, Runge-Kutta 
integration method (subroutine RUNGE). The number of integration steps has 
been chosen by making calculations of the penalized performance index and its 
derivatives for different step sizes and choosing the largest stepsize which 
gives reasonable accuracy. This has been done to minimize the computation time 
and to keep the integration as far away from round-off error as possible. 

The integration subroutine makes repeated calls to the subroutine DERIV 
which contains the differential equations of motion of the vehicle and the 
differential equations for the inequality constraints. Here, the tables 
containing Che controls are read; the aerodynamic, propulsion, and mass char¬ 
acteristics are determined; and the right-hand sides of the differencial 
equations are computed. The aerodynamic characteristics are obtained from a 
call to Che subroutine AERO which reads the aerodynamics tables. To obtain 
the atmospheric properties,>AER0 calls subroutine ATM62. 

Subroutine SLINl performs linear interpolation of a table of data points. 
It has been provided to read the control tables and the aerodynamics tables. 

The last subroutine supplied by the user is SGX. Here, the derivatives 
of the performance index and the constraints are computed by central 
differences. Hence, two function evaluations are made for each derivative 
computed. 
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SECTION VI 


SOLUTION OF THE OPTIMIZATION PROBLEMS 

The general procedure which has been followed in Che solution of the 
optimization problems Is to find a nonoptlmal, constant control trajectory 
which gives a reasonable value to the performance index and is somewhat near 
the desired final conditions. This is done by varying the parameters 
(final tine, angle of attack, bank angle, and ignition time if necessary) 
systematically until a reasonable trajectory is obtained. Then, these 
parameters are used as Che initial guess of the optimal trajectory. 

6.1 Reentry Problem 

For the reentry problem, a reasonable sec of nominal parameters has 
been found to be t^ - 3200 sec , a ■ 20 deg , and u ■ 60 deg . The 
results from integrating the equations of motion for these values of the 
parameters are presented in Table 5. Note that the nominal crossrange is 
•• 26.4 deg and that the final condition is given by 

« hj/100,000 - 1 - 0.136 (62) 

The main program contains all of the input quantities and ^s presented 
in Appendix C for the reentry problem (IPROB • 1). The iteration process 
starts from the nominal path (IRS ” 0), and no more than TMAX •• 1200 sec of 
computer time are to be used on the first run. To allow the program to reach 
the point where the time check Is made, the external tine limit has been set 
at 1500 sec. A total of eleven parameters (N • 11) is being used, of which 
one la the final cine, five (NA ■ 5) are the ordinates of the angle of attack 
table, and five (NM ■ 5) are the ordinates of the bank angle table. There is 
one constraint (M“ 1) which is an equality constraint (K ■ 1). Next, the first 
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Table 5. Nominal Path for Maximum Croasrange 
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parameter X(l) is set equal to the nominal final time, X(2) through X(6) 
are set equal to the nominal angle of attack history, and X(7) through X(ll) 
contain the nominal bank angle history. Also, the values of the normalized 
time at each of the nodes are the same for each table, TTA(I) and TTM(I) , 
and are given by 0.0 , 0.25 , 0.5 , 0.75 , and 1.0 . The unconstrained 

minimization will terminate when the relative change in each parameter between 
two iterations is less than EPS(I) - l.E-4 * X(I) . On the other hand, the 
constraint satisfaction iteration will terminate when the constraint residual 
satisfies the relation C(l) AKMIN * C{M + 1) where AKMIN - l.E-4 and 
where C(M + 1) is defined below. Since C(M + 1) “ 0.136 , an accuracy 
of l.E-5 is being requested In the constraints. 

The expected change in the performance index on the first iteration is 
DFN = 0.5 . MAXEN is set equal to a large number (10,000) because the program 
is now Interrupted on a time basis rather than a function evaluation basis. 

Every iteration in VFOIA is to be printed (IPRl » 1) as is every iteration 
in VA09A (IPR2 “ 1). The work space W(I) is made to have IW - 2500 words as 
suggested by Appendix A. The optimization code will calculate the weight 
(MODE = 1) using the constraint scale factor C(M + 1) = 0.136 . To modify 
to code for interruption, it is now necessary to input the constraint scale 
factors in C(M + I) rather than in C(I) as stated in Appendix A. Also, 
the parameter AK <• 1E60 must be set here as well. The remaining quantities 
and the rest of the main program direct the writing on tape and reading from 
tape if the computation is interrupted (CPU time > TMAX). 

A summary of the results is presented in Tables 6 though 9, and the 
characteristics of the converged trajectory are shown in Figures 1 through 3. 

Four cycles have been needed to achieve the desired accuracy in the final condi¬ 
tion, It is interesting to note that a good trajectory is obtained after one cycle 
and that the remaining cycles are used to satisfy the final condition better. The 
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optimal angle of attack seems to be at Che value for maximum llft-to-drag 
ratio, although this has not been verified. With regard to Che roll angle, 

Che vehicle starts out at a high roll angle, eo that it dives Into the atmos¬ 
phere, and Chen gradually reduces the roll angle until It Is heading northward 
with a very small roll angle. The vehicle achieves a crossrange of approxi¬ 
mately " 47 deg and experiences a maximum load factor of around 1.5 g's. 

It should be noted that the weight has not been Increased during 

the optimization process. This means Chat, after each cycle, sufficient 
progress has been made in converging to the constraints. On the other hand, 

3^ Is changed during each cycle, and the Lagrange multiplier = ‘^1^1 
appears to be approaching a limiting value. This multiplier and the optimal 
controls can be used to generate Initial values of the Lagrange multipliers 
needed for starting Che shooting method. 

The complete computation required 1640 sec of CPU time (Cyber 170/750) 
and a total of 3000 function evaluations. This amounts to 0,55 sec per function 
evaluation, or 12.6 sec per function/derivative evaluation (23 function evalu¬ 
ations). The integration has been performed with 250 Integration steps, and 
it may be possible Co reduce Che time by using forward differences. However, 
the convergence characteristics would deteriorate. 

6.2 Plane Change Problem 

For this problem, a good Initial guess of the final time, the angle of 
attack, the bank anple, and the ignition time have been found to be t^ ■ 2733 sec , 
(1-40 deg , u * 60 deg , and t^ • 1624 deg . The corresponding trajectory 
Is shown In Table 10, from which It is seen that the orbital Inclination Is 
1^ ■ 22.2 deg and Che constraint residuals are 
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Noalnal Path for Maxinum Plane Change 
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.863n03E-2 


: h^/608,000 - 1 - 

Cj = V^/23.500 - 1 - .882155E-2 

Cj = Yf • ■109798E-1 (63) 

= y? - -0 

S ^ ^8 * 

The last two constraints are the Inequality constraint on the angle of attack 
(a ^ 40 deg) and an inequality constraint on the bank angle (p ^0). In the 
solution of the minimization problem, the bank angle at the final point tried 
to go to a large negative value. The last inequality constraint has been 
Included to keep this from happening. 

The set-up for the plane change problem (IPROB - 2) Is similar to that 
of the reentry problem. An additional parameter, the Ignition time, Is 
present making a total of twelve (N ■ 12). Here, there are five constraints 
(M *5), of which three are equality constraints (k ■ 3). At first, the values of 
EPS had been chosen to be l.E-4 * X(I) • However, because of the extreme 
sensitivity of the problem to the ignition time, it has beer, necessary to use 
EPS(I) ■ l.E-5 * X(I) . Finally, the constraint scale factors C(M + I) for 
the altitude, the velocity, the flight path inclination, and the ineouallty 
constraints have been set equal to .009 , .009 , .011 , .001 , and 001 , 
respectively. The first three values are rounded constraint residuals,and the 
last two have been chosen to produce relatively small weights. 

A summary of the iteration process is presented in Tables 11 through 14, and 
the converged trajectory is illustrated In Figures 4 through 6. After four 
cycles, the constraints are satisfied to four significant figures. The 
converged angle of attack history stays near 40 deg (maximum C ) at the 
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Table 12. Maximum Plane Change; Cycle 2 Results 
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Figure 5. Maximum Plane Change: 6 , , and h versus 




(deg) 



Figure 6. Maximum Plane Change: V , » and i versus 










high altitudes ^md dicraases to 13 deg (naxlnmin C /C-) at the low altitudes. 

Irf 1/ 

On the other hand, the bank angle Is near 70 deg during the entry part of 
the trajectory and decreases CO 0 deg during Che exit part. The maximum plane 
change Is 1^ > 33.7 deg , and the load factor does not exceed 0.5 g's. 

The CPU time for the four cycles la 4175 sec on the CYBER 170/750 
computer. This amounts to approximately 0.5 sec per function evaluation and 
10 sec per derivative evaluation. 

For this problem, the Lagrange multipliers v do not seem to be approach¬ 
ing a limit. A possible reason is Chat the problem is so sensitive to the 
value of c. that additional accuracy is needed in the numerical integration 

D 

to achieve better convergence accuracy. Also, In generating numerical deri¬ 
vatives, the same perturbation size is used for each variable. More consistent 
results can be achieved by tailoring the perturbation size to each variable. 

The behavior of the bank angle of the final point can be explained by 
assuming that the earth Is not rotating. The Initial point r.f the trajectory 
is located on the equator, and the vehicle la trying to put the velocity vector 
on the orbital plane with the highest inclination. As long as the longitude, 

6 ,1s less than 180 deg, the bank angle for maximizing the inclination is 
positive. However, once the vehicle crosses the equator into the southern 
hemisphere, it must bank In the opposite direction (p<0) to obtain additional 
Inclination. By restricting the bank angle to positive values, the maximum 
inclination should be obtained when 6^ > 180 deg and p - 0 . The optimal 
trajectory presented here has these features, as can be seen from Table 14. 

On the other hand. It should be possible to increase the inclination by 
relaxing the inequality constraint to u ^ -90 deg . 
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SECTION VII 


DISCUSSION AND CONCLUSIONS 

The maximum crossrange and the maximum plane change problems have been 
formulated as parameter optimization or nonlinear progranmlng problems. An 
existing code for solving the nonlinear programming problem vlth the augmenCed- 
Lagranglan method has been selected and modified slightly for solving these 
problems. The derivatives needed for the optimization method have been computed 
numerically. As a consequence, the entire program can be viewed as a model 
for solving virtually any trajectory optimization problem. 

A difficulty In solving trajectory optimization problems lies In the 
model used for the problem. For example. In performing the one-dlmenslonal 
search. It Is possible to change the parameters so much that resulting trajec* 
tories go to extremely high altitudes or extremely low altitudes. This has 
caused an exponential overflow In the atmosphere subroutine. In another case, 
the resulting trajectories ended up in a vertical dive, causing the numerical 
Integration to blow up because rf the singularity In the heading angle equation. 
Initially, this difficulty had been eliminated by removing the singularity, 
which means Increasing the number of differential equations by one and Increasing 
the numerical complexity of the program (time becomes a dependent variable). 
Later, by additional modification of Che optimization code. It has become 
possible to monitor the flight path Inclination and terminate the one-dlmenslonal 
search when y gets too small. Another difficulty arose by trying to modify 
the drag polar so that the Inequality constraint on the angle of attack could 
be eliminated. The optimization process has found a way to use the change 
to its advantage and has produced an unrealistic trajectory. Finally, in che 
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plane change problea, the bank angle at the final pi nt wants to go to a 
large negative value. At this point. It Is not understood what aspect 
In the model makes It useful for this to happen. Hence, an Inequality 
constraint has been Imposed on the bank angle to keep it nonnegative. 

Once the modeling problems have been eliminated, the reentry problem 
Is relatively easy to solve. This Is not the case for the plane change 
problem, as can be seen from the number of Iterations and the computer 
time required. Possible reasons Include the number of constraints Involved, 
sensitivity of the problem to Ignition time, and numerical accuracy, either 
In the Integration or In the computation of derivatives. With regard to 
constraints, the plane change problem is easy to solve if the bank angle is the 
only control and the altitude is the only constraint at the final point. Adding 
the other final conditions reduces greatly the progress achieved on every 
iteration. Also, Including the angle of attack as a control reoulres the use 
of the inequality constraint on a . 

The particular method for computing the derivatives, central differ¬ 
ences, has been chosen because It is easy to modify the model. (If the 
derlvsti''es are computed from differential equations, any change In the 
model produces a change in the derivative equations.) Although It has not 
been done, it is possible to cut the computer time roughly in half by only 
integrating over that part of the trajectory which Is affected by a per¬ 
turbation. Also, the computer time can be reduced further by using forward 
differences, or one function evaluation per derivative calculation. This 
reduction In accuracy will affect the accuracy with which the constraint 
can be satisfied. 
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The Banner In which the optialzation probleBs have bean solved is 
sufficient CO denonscrace the ability of the optialzation code to handle 
Much problems. Questions concerning the number of nodes needed to obtain 
an accurate solution and Che placement of these nodes are still open, but 
they can be answered with a aoderate saount of computer time. 
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APPENDIX A 

DOCUMENTATION ASSOCIATED WITH 
THE OPTIMIZATION CODE 
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Harwell :»ubfoutin« Library 


VFQ] a/aD 


1. Purnose 

To find the ninieum of e function r(j(} of n verieblee x subject to both 
equsllty constreints c^(£} ■ 0 i ■ i,2,..Tilc end Inequelity"constraints 
^ 0 1 a k 4^ (k a 0 or k a m is allowed but k must be 4 n)a 

Derivatives of all func'. Ions with respect to x must be provided, both the 
vector (<3F/eXp dF/dx2,....,3r/dx„) and the matrix whose i^^ column is 
(dc|/dx|, dcj^/<)x2,«*»>,dc^/dx^) for i a These functions and 

derivatives must be specified in a user subroutine VPOlB (see section 4). 

An initial estimate of the solution must be provided which need not be 
feasible. The subroutine allows advantage to be taken of the possibility that 
some of the constraints are linear, and also of certain other types of infor¬ 
mation about the problem, if available. If all the constraints are linear, 
the use of VPOIA is not most efficient, and one of the LA or VE routines 
should be used. lYie method is a penalty function - Lagrangian method (see 
section 8) and VFOIA calls VA09A to carry out the associated unconstrained 
miniraisations. 

2. Argument List 

CALL VP0U(N,H,K,X,EPS,AKMIN,DPN,KAXFN,IPR1,IPR2,IW,MODE) 

N An ir/rEGCR set to the number of variables n (N ^2). 

M An INTEGER set to the total number of constraints m (M i^l). 

K An INTEGER sat to the total number of equality constraints k. 

X A .,laL array of N elements in which the initial estimate of the solution 

must be set. VPOlA returns the solution 9C in X. 

CPS A REAL array of N elements, in which the tolerances for the unconstrained 
minimizations must be set. EPS(I) should be set so that CPS(I)/X(I) % 
AKMIN, roughly speaking. 

AKNIN A REAL number in which the relative error tolerance required in the 
constraint residuals must be set. VFOIA will exit when max|lc].(x) I / 
scaling factor for c^j ^ AKMIN for the active constraints jlj (see 
section 7). 

PFN A REAL number in which the likely reduction in F(x) must be sot. This 
is done in the same way as for VA09A. - see the VA09A specifIcetlon 

sheet. 

MAXFN An INTEGER in which the maximum number of calls of VFOIB on any one 
unconstrained ainimizatlon must be set. ;)oughly speaking 2 or 3 times 
MAXFN calls of VFOIB are likely to be made altogether. 

IPRl An INTEGER controlling the frequency of printing from VFOIA. Printing 
occurs every IPRl iterations, except for details of increases to the 
which are always printed. No printing at all occurs (except for error 
diagnostics) if IPRl • 0. All printing controlled by IPRl is suitably 

annotated. VFOIA 1 
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IPR2 An INTECER controlllnQ th« frequency of printing from VAO'V.. it>ri2 
should be set as described in the specificaticm sheet. 

IW An INTCOCR qivlng the amount of storage available In COHMCN/VrOliywf.). 
Set to 2Vi00 unless wishing to change the restrictions (see Section 1), 

MODE An INTtCXR controlling the mode of operation of VfX>lA. If any positive 
definite estlnste Is available at the heaelan matrix of the penalty 
function, aet |h00c| ■ 2 or 3, otherwlaa aet if10Dc| (see VA09A 
specification sheet). If eatlsiates of the and paranieters are 
available (eee section 8) set HOOE < 0, otherwise set MODE >0. A 
normal setting for e one-off Jbb with no Information availabla Is 
ncoe - 1. 

3. named COMMON areas 


Certain named COMHOM areas suet be declared end set on entry to VFOlA. 


CO«0N/vroiE/C(l50> 


cof#ioN/vroiF/ec(2S, so) 


COMfON/vro 10/T( 150) 


CO»tfON/VroiI/G2P(325) 


Set scale factors (>0) for the constraints In 
C(1),C(2C(M). Choose the magnitude of 
these scale factors to give an IrdIcAtlon of the 
magnitude of the constraints evaluated about the Initial 
approximation jc* If any oonacralnts are 
violated by an amount greater in aiodulus than 
that which Is set, then the setting la Increased 
aceordinoly. These scale factors are transferred 
to C(M.l), C(M»2),.....,C(2M) by vyOlA. 

Set the derivatives of any linear constraints on 
entry rather than In VFOlB. This Is the moat 
efficient and the numbers are not disturbed. The 
manner of setting is described In section 4. 

If MODE < 0 is used, then set the parameters 
01,62, .... 0* in T(l), T(2),...,T(M) and the 
parameters in T(M*l),T(M»2),...,T(2M). 

The meaning of these peramaters may be found In 
the report TP552 - see section 8. 

If ImooeI • 2 or 3 set the estimated heaslan 
matrix of the penalty function in C2P(!>,..., 
G2P(N*(N.l)/2). The manner of setting Is that 
described in the specification sheet of VA09A under 
the heeding MODE. 


4. The use r subr outine VPOIB 

The user must declare a subroutine headed 

3UBR0(JTINE VF01B(N,M,X) 

REAL X(l) 

CCMMON/VFOIC/P 
COMMON/VFO10/0(50) 

COMMON/VP0lE/C{150) 

COMMON/VFOlF/OCi25,50) 

This subroutine must ta)ce x, the vector supplied in X(1),...,X(N) and set 
P(x) in f'l C.(x),».« In C(1),...,C(M)] (dP/0x.,...<)F/dXn)x In G(l),... 

C(N)( and aet Tic^/dX'j,,..tic^/JXn)x in GC(1,I),...,GC(N,I) for aTl I - 1,2,•• 


2 VFOIa 
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[t^xcepllni) the linear constrairits wl ii-ii shnita r.et on entry, as the 
numbers are constant]. Some time m ■ .ilso ho .-. .rf—i if required by also Inclu¬ 
ding COMM0N/vr01G/T(150) and by not -'valua-ii.4 GC<N,I) when 

C(I) ^ T(I) for any I > K. Note ttmi. the opi.lmm- values F(y), (dF/ds^,..., 
df/dxp, etc. are left in these name'l Common ,->r.'-.:. on exit from VFOIA. Note 
also that In the double precision version t.iio usci subroutine name is VFOIBD 
and there is a 0 appended to the named COMMON areas. 

5. Redefining named COMMON areas 

Local storage for VFOIA is through named COMMON areas. These have been 

sec on the assumption that k 4 and M ^ SO. If it is desirwd to r«nove 

either or both of these restrictions, then it is necessary to increase the 
storage available in some or all of these areas. This can be done by defining 
the named COMMON areas in the users MAIN with the increased storage settings, 
in which case the extra storage will be effective throughout the whole program. 
The complete list of named COMMON used by VFOlA and the corresponding values 
of N and H are as follows. 


COMMON/VFOIC/F,M,K,IS,MK, NU 
D/G(50) 

E/C(150) 

r/GC(25,50)* 

G/T(150) 

H/GP(50) 

I(G2P(32S) 

J/V(50) 

K/WW<150) 

L/W(2500)** 

M/ZZdOO) 

N/LT(100) 


independent of N and M 
2N 
3M 


N,M 

3M 

p (p • max(M,N» 
N*(N*l)/2 


?P 

?M 


Notes: ^ 

* On increasing M, when N < 25, redefine GC with bounds (N,M) net (25,M). 

VrOlA has been coded under this assumption, as it requires less storage. 

(VFOIA treats GC as a single suffix array). 

••For M very large, p^ storage locations may be prohlbltlvely large. However 
it is very unlikely that this amount of storage will actually be needed by 
VFOIA. (no problem has yet been encountered for which more than (?N)2 locations 
have been required). Hence in these circumstances, either declare w with 
(2N)^ locations (or whatever can be spared), and set the integer IW to this 
number in the calling sequence of VFOIA. If mere locations are required, 
then VFOIA will stop and print out the storage required. 

6. General 

Use of COMMON: 

Workspace! 


other routines; 


Inout/Output ; 


Restrictions: 

VFOU 3 


named COMMON only - see section 5 . 

5K words unless N or M is redefined, when it is not 
usually more than •'4 ']n2 ♦ nm » 0(max (M,N)) words. 

Cslls VFOlO (user subroutine), VFOlZ (private), VA09A, 
MCllA, MCllE. VA09A calls MCllB in addition. 

No input, all output on stream 6 (line printer), oon- 
trolied by user througn IPRl and 2PR2. 

N ^ 25 and M ^ 50 but can be lifted - see section 5. 


53 


Prtte of routinaj 


August 1973 


7. Accuracy 

It may be that VFDIA Is unable to achieve the accuracy requested In the 
parameter AKMIN. In this case a diagnostic is printed. To find the cause 
of this, first examine the print out of the VA09A Iteration. If this is 
anomalous (gd 0 for Instance) suspect a mistake In the programming of VPOli) 
particularly in obtaining derivatives. If VA09A seems O.K., then other 
possible causes are (i) there is no feasible point (in which case -»•> and 
Ci a const d 0), (ii) EPS has been set too large relative to AKMIN, (iii) 

AKHIN has been set to saiell relative to the machine precision, (iv) the 
problem is too lll«>cendltlonad. 

8. Method 

That described in R. Fletcher. "An Ideal penalty function for constrained 
optiinlzation", C.S.S.2. December, 1973. The penalty function for 
inequalities is 


d(x,j5,ff) ■ F(xi r ^ 

and each iteration Involves minimizing ^ix) for fixed After each itera¬ 

tion the 8 and £ parameters are varied so that the sequence of minima 
tends to The soTution of constrained problem. The value of the producT^ 

0^0^ tends to the i^^ Lagrange multiplier of the problem. Any information 
about Lagrange multipliers, or about the hessian of d can usefully bo 
incorporated. 

Convergence Is guaranteed (in exact arithmetic) and this implementation 
of the method can be expected to converge at a second order rate. 


October, 1973 


VFOIA 
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Harwell Subroutine Library 


VA09A/AD 


1 . Purgose 

To find the mlnlnjin oT a flinctlon F(x) or several variables, glveit 
that the gradient vector (3F/0Xj,OF/&x2,.T.,OF/^Xj^) can be calculated. 

The subroutine replaces VAOtA to which it is superior in various v«ys (see 
section S), and should be used whenever derivatives can be evaluated readily. 
It should however not be used either if storage space is at a premium (use 
VA08A) or if the function is a sum of squares (use VA07A). TT»e subroutine' 
complements VA06A, the latter requires four times the storage, and some 
comparisons (R. Fletcher, A.E.R.E. Report, R712S (1972)) indicate that 
VA06A is marginally slower and more affected by round off error. As VAOM 
Is more difficult to use, it is suggested that VA09A should be 
used in the first instance on af\y problan. If VA09A fails then VAO^ should 
be tried as it is guaranteed to converge if the effect of rounding errors can 
be neglected. 

2. Argument List 

CALL VAOaV(FUNCT,N,X,F,G, H,W, DFK, EPS,MODE,IKXFN,IPRINT,lEXIT) 

FUNCT An IDENTIFIER of the users subroutine - see section 3. 

N An INTECaER to be set to the number of variables (N 5. 2). 

X A REAL ARMY of N elements in which the current estimate of the 

solution is stored. An initial approximation must be set in X 
on entry to VAOM and the best estimate obtained will be 
returned on exit. 

F A REAL number in which the best value of F(x) corresponding to 

X above will be returnee. 

G A REAL ARRAY of N elonents in which the gradient vector 

corresponding to X above will be returned. Not to be set on 
entry. 

H A REAL ARRAY Of N*(N+i)/2 elements in which an estimate of the 

2 

hessian matrix 9 F/(3x.ox ) is stored. The matrix is represented 

ij J 

in the product form LDL where L is a lower triangular matrix with 
unit diagonals and D is a diagonal matrix. The lower triangle 
of L Is stored by columns in H excepting that the unit diagonal 
elements are replaced by the corresponding elements of D, The 
setting of H on entry is controlled by the parameter M)DE (q.v.). 

W A REAL ARRAY Of 3«N elements used as working space. 


VA09A 1 
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Dm A REAL number which nuii. I>t' --ct so as to g.LVC VA09(L an 

estimate of the likely rcuui cion to be ol>t:tined in F(x). 

DFW is used only on the fii't iteration so an order or 
magnitude estimate will ^lUll'ice. The inrormetion can be 
provided in different ways depending upon Uic sign of DFN which 
should be set in one of the followlrig v«ys: 

Dm>0 the setting of DFN itself will he taken as 
the likely reduction to be obtained In F(x). 

DFMsO it will be assumed that an estimate of the 

minimum value of P(j{) haa bean set in argument 
F, and the likely reduction In F(x) will be 
computed according to the inltial~functlon value. 

Oim<0 a oultlple |Dm| of the modulus of the initial 
function value will be taken as an estimate of 
the likely reduction. 

EPS A REAL ARRAY of N elements to be set on entry to the accuracy 

required in each elanent of X. 

M)DE An nfTEGBl which controls the setting of the initial estiamte 
of the hessian matrix in Che parameter H. The following 
settings of KCDE are permitted. 

MOCE^I An estimate corresponding to a unit matrix 
is set in H by VA09A. 

M0I£=2 VA09A assumes that the hessian matrix itself 
has been set in H by columns of its lower 

T 

triangle, and the conversion to LOL form 
is carried out by VAO^. 1he hessian matrix 
nust be positive definite. 

MDCE^ VA09A assumes that the hessian matrix has been 

set in H in prodict form. This is convenient when 
using the H matrix from one problem as an Initial 
estimate for another, in vJilch case the contents 
of H are passed on unchanged. 

MVXFN An INTEGER set to the maxiinum number of calls of FUNCT permitted. 

IPRlffT An IVTEdR controlling printing. Printing occurs every 
|IPRINT| iterations and also on exit, in the form 

Iteration No, No of calls or FtJNCT,lEXIT (on exit only) 

Function value 

X(1),X(2),...,X(N) 8 to a line (5 in VAOBAD) 

G(1),G(2),...,C(N) 8 to a line (5 in VA09AD) 

The values of X and G can be suppressed on Intcniediate 
iterations by setting lPRiNT<0. All intermediate printing can 
be suppressed by setting lPRINT=MAXFFH’*i. All printing can be 
suppressed tiy setting IPRINT=0, 


Z yA09A 
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lEXIT An INTEGER giving the p -ion for exit from VAOtti. Thle 
will be aet by VA09A aa ■ollowe 

lEXIT^ (M0DE=2 only). The eatunate of the heaaian 
matrix ia not positive definite. 

IEXIT=I The normal exit in w ich |DX{1)|<EPS(I) for 
ali 1=1,2,...N, where DX(I) is the change in 
\ on an iteration. 

T 

IEXI'ha2 G DK^. Not possible without rounding error. 
Prebable cause is that EPS is set too fliall 
for computer word length. 

IEXIT=3 FUNCT called kKXFN times. 

3. User Subroutine 


The user nust provide a subroutine headed 
SUBRCMTINE XXX(N,X,F,G) 

REAL X(1),G(I) (REAL*8 in VAOSAD) 


v4iere XXX is an iaentifier chosen by the user. 

This subroutine should use the variables x supplied in X(l), 
X(2),...,X(N) to evaluate the function and gradient vector and place 
them in F and G(1),G(2),...,G(N) respectively. . XXX mist be passed to 
VA09A as VAOfiA'a first argument, see section 2, and appear in an 
EXTBINAL stateatent in the program that calls VA09A>. 

4. General 

Use of OOMGN ; none 

Workspace: N*(N+l)/2 words * 4N words provided by the user in 

H and W. 


Other routines; calls MC11A, B, E. 

Input/TXitput: controlled by tho user through IPRINT, All 

output is on stroam B (line printer). 

Restrictions; none 


Systan dependence; none 

Date of routine: April, 1672. 
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5. Method 


The method used is a ouasi>^le«rton method described by Kletcher 
(Ooaqxjter Journal, Vol. 13, P,3I7, 1870), and is a modirication of 
earlier methods of this tyr such as that implemented by VAOIA, The 
method is superior to that • < VAOIA on three counts* 

(1) It uses a formula to update the hessian approximation H 
which has proved to be more efficient and reliable. 

(2) It uses a 'crude' line search %vhich has been shown to be 
more efficient than an 'accurate' line search. 

T 

(3) It represents H by the product LOL , which enables the 
positive definiteness of H to be guaranteed, even in the 
presence of round-off error. 





HirMell Subroutine Llbrery 


MCI I A/AD 



i 

I 


t< Purpoee 

MC!IA la a subroutine for use In probleas which involve the addition 
or subtraction of rank one aatrices <r xt* to positive definite-or seal- 
definite synnetric aatrices A stored In factored fom A s LOL^, such that 
the resulting N x N natrix 

*** T 

A =‘A * <T Xt 


is also known to be positive definite or seal-definite. Note that L is 
lower triangular with *hd D is diagonal with d^ ^0. Apart froa 

its obvious use in updeti.ng satrices which reaain strictly positive 
definite, yCilA can be used (i) to acauailate a sua of rank one teias into 
an initial aatrix Mi, (li) to carry out projection end allied 
opcratione on A which reaice the rank, and (lii) to update aatricee A of 
rank k < n where it ie known frcn other oonelderations that the rank 
rran) ns unchsnged. All these operations preserve the correct rank and are 
not seriotuly affected by round-off error. the aethod is that deecrlbad 
by M.J.D. Powell and R. Fletcher (|975), 'On the updating of LOL* 
factoriseations', T.P. 510. 

the aecrix A is rq>resented using the ailnlaal storage of N*(N-fl)/2 
eiraents where N la the diaension of the problea. Tb facilitate 
operating with A, a nuaber of independent eubroutlnes have been provided, 
written as entry points MClIB/C/.../F. These opereUons include reducing 
e aatrix to its factors, ojiupiying out the factors, operating with the 
rectors of A on a vector j to obtain either As or A~'z, and replacing the 
faccora of A by the aatrix A~'. These fscillties are described in aore 
detail in section 4. 

2. Arffissfp^ Upt 

CALL IC11ACA,N,Z,S10,W,1R,M(,EFS) 

A A REAL one disensionai array of N*(N'fl)/2 eienents in vhich the aatrix 
T 

AsLOL liiBt be glvfln in factored fons. Ihe order in vdiici) elonenti of 

L and D are stored is <»,.«2t'Sl»***'*NI'‘‘2'^32"‘”*N2'. 

^1'^ ^ factors of the aatrix 

T 

A E A az j will overwrite those of A on exit. 

N An iNIEtZR (N>l) which aust be set to the diaenaion of the problea. 

Z A REAL one diaenaional array of N eieaents in which the vector z auet 

be eat. Tha array Z Is ovirwrittan by the routina. 

810 A REAL variable in which the scalar <t lUst be set. SIO is not 
restricted to i1., but if SIC<0 then it sust be known froa other 
considerauons that X l# positive definite or seal-definite, apart 
frca the effects of round-off error. 
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W A REAL array of N eiententa. If SIG>0 then W is not used, and the 
nsBM of ai\y one dimensional array can be inserted in the calling 
sequence. If SIC<0 then W is used as work space. In addition for 
SIG<0 it stay be possible to save tijae by setting in W the vector 
^ defined by Lvb^* Ihe ways in which tMs can occur are described 
under MC belowT 


IR An INIECOl to be set so that |1R| is the rank of A. If the rank of A 
is expected to be different from that of A, set IR 4 O. On exit from 
MCIIA, IR(d^O) will contain the raak of A. 


MC An INTCCXR to be set only when 8IG<0, as foliowa. If the vector 
V defined by has not been calculated previously^ set MCaO. 

7r MCI IE ibs been used previously to calculate A'^z, then v is a by¬ 
product or thic calculation and is stored in the W~paraBeter of MCliE. 
In tias case transfer v to the W parameter of MCIIA and set MCsf. 

If z tub been calculated as z » Au for some arbitrary vector u using 
MCiTd, then again v is a by-prodict of the calailation and is~ 
available in the W^parameter of MCI ID. In this case (or spy other 
in which v is known) set v in the W parameter of MCIIA and set MCs2. 

EPS A REAL variable to be set only when SIG<0 and A is expected to liave 
the same rank as A. ^ In ££ct^in ill-conditioited cases a noi>-zero 
diagonal elesmnt of D (XsLOl/) might become so smell as to be 
indeterminate, 1V«o courses of action are possible. One is to 
introduce a small perturbation in order that X keeps ttie same rank 
as A. Ihis is the normal course of action aivj is achieved by setting 
EPS equal to the relative machine precision e. e is I0~^ in single 
length arithstetic and ^ |0~*^ in double length on the IBM 370, The 
other course of action is to let the rank of 7( be one less than the 
rank of A. Thia is achieved by setting EPS eiiual zero. 

3. Oeneral 

y2£^£_OOSSON: noiw 


Workspace : 


N*(N'fl)/2 ♦ N -f N words provided by the user in A,2 end W. 
If SIG>0 W f>eed not be supplied. 


Entry points ; 


MClIB/C/.../F - see section 4, 


Other routines; none. 


InoutAXitnut : 

Tim^: 


none. 

2 

One cell of MCIIA requires n multlpllcdtions, unless 

2 

SI0<0 and M(aO when the figure is l^i n . One call of 

3 / 

any of MCI IB,C or F requires n /6 suitiplicatlcns. One 

2 

call of either MCI ID or E requiree n /2 sultipiications. 


Reetrlctions ; none. 

Sys tem dni>t‘ndunce ; none. 

Date of rouane: Jaiuery, 


1 

•1 

1973. ' 

.1 







4, Outer Entry polnf 


1 


Other entry points are provided to facilitate operating tviUt A 
v«hlch IS stored in compact form. In all of these A is a R£AL one 
dimensional array of N*(N-fi)/^ clementSi and is an integer set to the 
dimension of the problan. Each entry point la essentially an independent 
si'iroutine, and could be taiten out and written as such if desired* 

NCI IB - factorize a positive definite sysaaetric matrix given in A. 

ONLL MCnB(A,N,IR) 

A Mist contain the elesients of A in the order a,,«a, ,a .,a 

•n!.V..*-.-"*,” 

successive columns of its lower triangle). On exit A will be over¬ 
written by the factors L and D in the form described in section 2 under A. 

N Order of the matrix A. 

IR An IN'IECER set by MCtiB to the rank of the factorization. If 

the factorization has been performed successfully IftsM will be set. 

If IR<N then the factorization has failed because A is not positive 
definite (possibly cbe to round-off error). In this case the 
fsctors of a positive semi-definite matrix of rank IR will be found in 
A. However the results of this calculation are unpredictable, and 
M211B should not be used in an attempt to factorize positive eaai- 
definite matrices. 

T 

MCHC - imiltlply out the fsctors LDL to obtain A, 

CALL MC1tC(A,N) 

A Mist be sec in the factored form described in section 2 under A. 

On exit the factors will be overwritten by the explicit mstrlx A, 
the order of the dements being that described for input to NCI IB. 

N Order of the imtrix A. 

MCUD - calculate the vector z*=Az y¥here A is in factored form, 

CALL MCIID(A,N,Z,W) 

A Mist be set in factored form. 

N Order of the matrix A. 

Z A REAL array of N elements to be set to the vector z. On exit 
z contains the vector ” 

W A REAL array of N elements which is set by MCI 10 to the vector v 
defined by Lv=^. If this vector is not of interest, replace 1? by 
Z in the calling sequence to obviate the need to supply extra storage. 
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tCl IK - calculate the vector £*=sA”** wliere A la in factored fom. 

GALL IC11B(A,N,Z,N.IR) 

A Miet be set in factored fora. 

N Order of the natrix A. 

Z A REAL array of N elesienta to be set to the vector g. On exit Z 
contains the vector z*sA~^g, 

W A REAL array of N elements wliich is set by (CUE to be vector v 

defined by L;^^. If this vector is not of intereati replace N by 2 
in the calling sequence to obviate the need to supply extra storage. 

IR An INIEOER which Bist be set to the rank of A, 

iClIF - calculates the e)q>lici t matrix A~* froai the factors of A. 

CAU IC11F(A,N.1R) 

A Uist be set in factored form. ICMP will overwrite this by the 

elesients of the inverse matrix A~’. in the order ‘rt**2i***"SN ** 
for MCI IB. 

N Order of the matrix A. 

IR An INTEGER which must be set to cne rank of A. 

Notes. (i) ICMF should not be used to solve equations, in stilch case 
MCME should be used. ICnF is intended for applications in which the 
explicit elements of A"^ oust be examined, for exasq;>le in the use of 
variance>covariance matrices. (ii) MCME and F will RETUIM without 
calculation unless IRsN* 


m 

i 


January, 1973 
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Harwell Subroutine Library 


VE04A/AD 


1, Purpoae 


To find the value of ^ 


*1 


wtiich nunlnilzea a '-.jadratlc funcUon 


QC]() of n variables subject to upper and lower bounds on aooe or all of 
the variables. QCj) is defined by 

n n n 

Q(«) ■ - Ij . Y' - ^•*l’‘l 

1=^1 j = l 1 = 1 


where A la synnetrlc , and the bounjs are of the form I^ax^au^, It 

la permissible to let If required, and A is not restricted to being 

positive definite. The subroutine calculati^s the solution js = £, the mininuffl 
value Qig), and the gradient g(^) (note 4 (x)=Ax-Ji), This problem is u 
special case of quadratic programming for which a subroutine VE02A exists. 
VE04A is more efficient and tellable for sol' ing problems of this special form. 

An applicatioii of the subroutira! Is to (weighted) llnrar least squares 
dsts fltUng subject to bourds. If it is required to minimize 


S(2c) = (Bi5-y)‘''w(B4-y) 


( 2 ) 


subject CO the above bounds, where b Is an mxn matrix m>n and W Is an mxin 

diagonal imirlx of weights tnen set A=2B^WB arfl ^ = 2B^»yy. Statistical 

calculations for this problem are described in section 5, including an additional 

entry point to enable the variance-covariance natrlx to be calculated. 

2. The Arioinent .list 
♦ 

CALL VE04A (N,A, IA,B,BL,OU,X,(J,Lr,K,C) 

N an irffCCER which must be set by the user to the nuAber of variables. 

A a REAL, two dlneiw lonii 1 urrqy, iNidi diiiL-nslon atlroBt N; Che 

elements in the upper tciansl'' A(l,J) TaJcn must be set by the 

V30U 1 
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user to tite correipoivllng In (I), ond «dll rcnBin untouched 
the eubiDuUne. Elenents A(1,J) N>I>J arc used as Morking space. 

lA an INTEGER giving the first dinenaion of A in the stateiKnt which 
assigns space to A. 

B a REAL array of at least N eicncnts. The user oust set B(I) to the 
in (I). B is not oveivrlttsn by VE04A. 

BL a REAL array of at least N elenents. The user aaist set BL(1) to 

the lover bound on the variable. If the bound is non*existent, 
set to a very sisali nsiber like -1E7S. BL is not oveiwritten by 
VE04A. 

BU a REAL array of at least N elenents. The user must set Bu(I) to the 
upper bound on the I^ variable. If the bound is non-existent, 
set u^ to a very large number. BU is not overwritten by VE04A. 

X a REAL array of at least N elenents. VE04A returns the suiutlon 
InX(l). 

Q a REAL variable in which VECHA returns (he value of Q(^). 

LT an INTEGER array of at lesat N eleaents, aet by VE04A to a 
permutation of the integers 1,2,...,N (see K and C below) 

K an INTEGER set by VE04A to Uie nusber of free variables at the 
solution j (those not on Ueir bounds). These are the variables 
LT(l), LT(2),...,LT(K). 

C A REAL array of at least 3*N eleaents, C(1),....,C(N) are set by 
VE04A to the gradient 4((}., C is indirectly addressed so that 
C(I) contslre the gradient with respect to the LT(I) variable, 
whence G(1 C(K) will be found to be zero. CCN'fl 

CON) are used by VE04A as working ^sce. 


2 VdOLA 


64 





3. 


Goirrul Infonmtlon 


Use of COkMON ; none 

2 

Workspace : approx Jjn wrds (fair of A) + 2n words In G and n iiitegera 
in LT. 

I her routines ; none 
Entry points ; VE04B-8ee section 5 
Input/Output : none 
System dependence : none 

Timing. The time required depends upon how many free variables k there 
are at the solution. Typical f^ ures for (k/n, nunher of 
multiplications) are (. 1 , n^/|2), (.S^nVe), (.75,nVlf). 
Restrictions : none 
Date of Routine : April 1973, 

4, Mh thod 


That of Fletcher R, and Jackson, M.l*., (1973), "Minimization of a quadratic 

function of many variables subject only to lower and upper bounds", T.P.528, 

This method contilnes generality (any A), efficiency (times comparable to those 

T 

required for a factorization of A) and stability (uses partial LDL factorizations). 
5. Statistical Calculations 

When a sum of squares is being mininitzed as in (2), then certain statistical 

quantities can readily be calculatr-i Firstly, of course, S(§) is given by 
T 

Q(^) . _y h[y. If it is assumed that tit bound variables located by VE04A are 
exact in the underlying model, thc(\ an,estimate of the residtal variance is 
given by S(^)/(m-k). To estimate vanatscs and covariances, an additional 
entry point is provided. This calculates {aJ * »>1iere lAl Indicates the submatrix 
iA^ji where i and J index only the free variables. Tic appropriate variance- 
covariance matrix for the free variables is then {a!" '. Estimates of 

stanilanl deviations of the free variables are given by the square roots of the 
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diagonal clinvntH of this matrix. Because the bound variaoies are known 
exactly. Uvy have aero variance and covariance. 

The entry point VE046 is esaentisily written aa a separate subroutine. 
It calculates and is used as follows: 

CALLVE04B (N,A,U,G,K) 

N,A,IA,G and K aaist be paased on unchanged after exit from VE04A. 
VE04B sets the following: 

A The off-diagonal elesients of (Aj~* are set in A(I,J) for 

J<1«K. The elements are indirectly addressed so that A(I,J) 

contains {a}"* where r=LT(l) and saLT(J). 
r .a 

C The diagonal elesmnts of iAi~' are set InG(N^l).. C(N-t-K). 

They are again indirectly addressed so thatG(N-^I) contains 
-1 

where r=LT(I). 


2nd May 1973 
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VF01A PAGE 1 


SUBROUTINE VFO1A(N.N.K.X.EPS.AKHIN.OFN.NAXFN.IPR1.IPR2.lU.NODE) 
REAL XCD.EPSd) 

COMMON/VFO1C/F,MM,KL.IS,KK,MU 
COMM0N/VF01D/G(50) 

COMMON/VF01E/C(150) 

COMM0N/VF01F/GC(1250) 

COMMON/VFO1G/T(150) 

COMMON/VF01H/GP{50) 

COMMON/VFO1I/C2P(325) 

COMMON/VF01J/V(50) 

COMMON/VFO1K/WW(150) 

COMMON/VFO1L/W(2500) 

COMMON/VFO1M/ZZ(100) 

COMMON/VFOIN/LTdOO) 

COMMON/TIME/TMAX,TO,T1,IRS,MINS,AK,ITM,IFN 
COMMON/INON/INOM 
EXTERNAL VF01Z 
DATA BOUND/1. E-t-S/ 

1000 FORMATOOI**) 

1001 FORMAT(8E15.7) 

MUsMAX0(25,N) 

IF(M.GT.50)NUsN 

IXsN 

ICSsM 

ICBsM^M 

ISsM 

IL=IS+M 

IPsM 

ILT=M 

NN=N»(N+1)/2 

MMsM 

KL=K 

R=1. 

MK=0 

INOMsI 

CALL VF01B(N,M.X) 

DFsDFN 

IF(DFN.LT.OEO)DF=ABS(DFN»F) 

IF(DFN.EQ.OEO)DFsF 

IF(DF.LE.0.>DF=1. 

IF(HODE.LT.O)GOT05 
DO 2 1=1,M 
CC=C(I) 

IF(I.GT.K)CC=AMIN1(CC,0E0) 

IF( ABS( CC) .GT.C( ICS4.I) )C( ICSt.I) = ABS(CC) 

2 CONTINUE 
DO 3 Isl.M 

T(ISi.I)=2E0»DF/C(ICS*I)**2 

3 T(I)=0. 

5 CONTINUE 
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IF(IPR1.EQ.O)GOTOi< 

IF( MOD( HINS. IPR1). NE.0 )GOTOil 
PRINT 1002 

1002 FORNATCOENTRY TO VF01A'///*OCONSTRAINT SCALE PARAMETERS ARE') 
PRINT 1001,(C(ICS»I),1=1,M) 

U CONTINUE 
HD=IABS(HODE) 

8 CONTINUE 

DO 9 Isl.NN 

9 W(I)=G2P(I) 

IF(IPR1.EQ.0)G0T07 
IF(M0D(MINS,IPR1).NE.0)G0TO7 
PRINT 1003,MINS 

1003 FORMAT(////'OOUTER ITERATION NUMBER IS'.I3> 

PRINT 10011 

100«» FORMATCOXCI)') 

PRINT 1001,(X(I),1=1,N) 

PRINT 1005 

1005 FORMATCOTHETAd)') 

PRINT 1001,(TCI),1=1,M) 

PRINT 1006 

1006 FORMATCOSIGMAd)') 

PRINT 1001,(TdS+I),1=1,M) 

7 CONTINUE 

IFdRS.EQ.1) GO TO 82 
ITN=0 
IFN=0 
82 IRS=0 

CAa VA09A(VF01Z,N,X,PHI,GP,W.WW,DF,EPS,MD,MAXFN,IPR2,1EXIT) 

IF(TI-TO.LT.TMAX) GO TO 80 

DFN=DF 

DO 81 1=1,NN 
8i G2Pd)=Wd) 

RETURN 
80 CONTINUE 
INOHrl 

CALL VF01B(N.M,X) 

MD=3 

AKK=0. 

DO 10 1=1,H 

IFCI.GT.K.AND.Cd) .EQ.O.) Cd) = l.E-8 
CC=Cd) 

IFd.GT.K.AND.CCI).GE.TCI))CC=AMIN 1 (CC,OEO) 

Td) = Td)»T< IS^I) 

WWd)=ABS(CC)/CdCS+I) 

IF(WWd) .GT.AKK)AKK=WWd) 

10 CONTINUE 
IFdPR1.EQ.O)GOTOl6 
IF(M0D(MINS.IPR1).NE.O)GOT016 
PRINT 1007 
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1007 PORMATCOEXIT FROM VA09A>/'0LAGRANGE MULTIPLIERS USED IN VA09A'} 
PRINT 1001,(T(I),Is1,M) 

PRINT 1008 

1008 FORMATCOLARGEST SCALED CONSTRAINT VIOLATION'/ 

1' THIS ITERATION, BEST ITERATION') 

PRINT lOOl.AKK.AK 
PRINT 1009 

1009 FORHATCOCONSTRAINT RESIDUALS') 

PRINT 1001,(C(I),1=1,M) 

PRINT 1010 

1010 FORMATCOSCALED CONSTRAINT VIOLATIONS') 

PHINT 1001,(WW(I),1=1,M) 

16 CONTINUE 

IF(IEXIT.EO.0.OR.IEXIT.EQ.3)GOTO20 
IF(AKK.LE.AKMIN)GOT020 
IF(AKK.GE.AK)G0T011 
DO 15 1=1.NN 
15 G2P(I)=W(I) 

DO 17 1=1,M 

IFd.GT.K.AND.Td) ,EQ.OEO.AND.C(I).GT.OEO)GOT017 
2ZdP+I)=-TdSfI)»Cd) 

IFd.GT.K.AND.2ZdP-*.I).LT.-Td))2ZdPd)=-Td) 

17 CONTINUE 
IF(MIMS,EQ.1)G0TO40 
GOTO18 

11 CONTINUE 

DO 14 1=1,H 

IF{ WWd). LE. AK. OR. CdCB«-1).GE. 4E0»V/Wd) )G0T014 
DS=9EO*T<IS*-I) - 
TdS>I) = lE1*T( IS+I) 

IFdPR1.NE.O)PRINT 1011,1,TdS+I) 

1011 FORMATCOSIGMAC ,13,') INCREASED TO ’,E15.7) 

DO 12 J=1,N 

12 V(J)=CC(d-1)»NUW) 

CALL MC11A(G2P,N.V,DS,V,N,N,DS) 

14 CONTINUE 

18 CONTINUE 

DO 13 1=1,N 

IF( ABS(Xd)-CdX-*-I)) .GT.EPSd) )G0T021 

13 CONTINUE 
PRINT 1013 

1013 FORMATCOREQUESTED ACCURACY NOT OBTAINED') 

20 CONTINUE 

DO 2012 1=1,NN 
2012 G2Pd)=Wd) 

IFdEXIT,EQ.O)PRINT 2000 

2000 FORMATCOMATRIX SET IN G2P BY USER IS NOT POSITIVE DEFINITE') 
IFdPR1,E0.0)RETURN 
PRINT 1012 

1012 FORMATCOBEST SOLUTION OBTAINED'/'--.,(0(1) ,1=1 .N) •) 
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PRINT 1001.F,<G(I).I«1.H> 

RETURN 
21 CONTINUE 

IF(AKK.LT.AK)G0T040 
DO 32 1=1,M 
32 V(I).T(IU1) 

COTO70 

40 CONTINUE 
HK>0 
KKsO 

DO 41 1=1,M 
T(IUI)=T(I) 

C(ICB4>I)=WW(I) 

IF(I.CT.K.AND.T(IUI).EQ.OEO.AMD.C(I).GT.OEO)GOT041 

KKsKK4>1 

LT(ILT+KK)=I 

GP(KK)s-BOUND 

IF(I.CT.K)GP(KK)=-T(IUI) 

VCKKIsBOUND 

ZZ(KK)s-C(I) 

41 CONTINUE 
IF(KK.EQ.0)GOT020 
DO 42 1=1 ,N 

42 CdX+DsXd) 

KKKaKK«(KK^1)/2 

IIaHAX0(KKK4-NN,KK>KK) 

IF(II.LE.IW)COTO50 
PRINT 2001,11 

2001 FORKATdOINCREASE STORAGE IN COMMON/VF01L TO',17,’ELEMENTS') 
RETURN 

50 CONTINUE 
IIsIW-KKK 

DO 53 1=1 ,ICK 
LI=LTdLT>I) 

DO 51 JJ=1,N 

51 X(JJ)=GC((LI-1)»NU-fJJ) 

CALL MC11E(W.N,X,DUM1.X.N.IDUH.DUH2} 

DO 53 J«1,I 
LJ.LTdLT^J) 

ZsO. 

DO 52 JJs1,N 

52 Z«Z4.X(JJ)»CC((LJ-1)«NU4.JJ) 

IIsII^I 

53 WdD.Z 
JJsIW.KKK 
11=0 

DO 56 1=1.KK 
DO 55 J=1,I 
JJiJJ^I 

55 WdI*J)»W(Jv') 
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56 IIsII+KK 

CALL VE04A(KK,W,KK,ZZ,GP,V.T,Z.LT,JJ.WW> 

IF(IPR1.EQ.C)G0T059 

IF(MOD(MINS,iPRI).NE.0)G0T059 

PRINT 1020,KK 

1020 FORMATdA,' ACTIVE CONSTRAINTS, NUMBERED*) 

PRINT 1000.(LT(ILT*I).1=1,KK) 

PRINT 1021 

1021 FORMATCOLAGRANGE MULTIPLIER CORRECTIONS FOR ACTIVE CONSTRAINTS*) 
PRINT 1001,(T(I).1=1.KK) 

59 CONTINUE 

DO 60 1 = 1 ,M 

60 V(I)=T(IL+I) 

DO 62 1=1,KK 
LIsLTdLT+I) 

V(LI)=V(LI)-fTd) 

Z=4E0*ABS( (Td)-2Z( IP+LI) )/Z2( IP>LI)) 

IF(Z.LE.1E0)G0T062 

DS=(Z-1E0)«TdS+LI) 

T( IS+LI) = Z*TdS>LI) 

IF(IPR1.NE.0)PRINT 1011.LI.TCIS+LI) 

DO 61 J=1 ,N 

61 GP(J)sGC((LI-1)*NUW) 

CALL MC11A(G2P,N,GP,DS,GP,N,N,DS) 

62 CONTINUE 
AKsAKK 

70 CONTINUE 

DO 71 1=1.M 

71 Td)sVd)/T(IS4.I) 

DO 72 1=1,N 

72 Xd)=GdX^I) 

DFs1E50 

HINS=MINS>1 

G0T08 

END 
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SUBROUTINE VA09A(FUNCT.H.X.F.C.H.W.DFN.EPS.MODE.HAXFN.IPRINT, 
1 lEXIT) 

REAL X(1).G(1).H(1),W(1).EPS(1) 

COMMON/TIME/TMAX.TO.T1.IRS.MIMS,AK.ITM,IFM 
COMMON/IGAMMA/IGAMHA 
IF(1PRINT.NE.0)PRIMT 1000 
1000 FORMAK//'ENTRY TO VA09A'/) 

NN=N«(N4.1)/2 

IG=N 

IGGsN-t-N 

ISsIGG 

lEXITsO 

IRsN 

DF=DFN 

IF(M0DE.EQ.3)G0T015 
IF(MODE.EQ.2)GOT010 
IJ=NN+1 
DO 5 1=1,N 
DO 6 J=1,I 
IJ=IJ-1 
6 H(IJ)=0. 

5 H(IJ)=1. 

GOTO15 
10 CONTINUE 

CALL MC11B(H,K,D1,D2,D3,IH.ID1,D4) 

IF(IR.LT.N)RETURN 
15 CONTINUE 
Z=F 

CALL FUNCT(N,X,F,G) 

IFN=IFN*1 

IF(DFN.EQ.O.)DF=F-Z 
IF< OFM. LT .0.) 0F=ABS< DF» F) 

IF(DF.LE.0.)0Fs1. 

20 CONTINUE 
DFNsDF 

CALL SECOND(TI) 

IF(T1-T0.GE.TMAX)G0 TO 90 
IF(IPRINT.£Q.0)G0T021 
IF<MOD(ITN,IPaiNT>.N£.0)G0TO21 
PRINT lOOl.ITN.IFM 

1001 FORMATCZ^IS) 

PRINT 1002,F 

1002 FORMAT((8£15.7)) 

IF{1PRINT.lt. 0/G0T021 
PRINT 1002,(X(l),1=1,N) 

PRINT 1002,(G(I),I»1,«) 

21 CONTINUE 
ITNsITH^I 
DO 22 1»1 .N 

22 W(IC>I’;»C<I) 
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CALL MClie(H,M,C.0^,W,IR,ID1,D2) 
GSsO. 

DO 29 TsI.N 
W(1S+I)=-G(I> 

29 GSrGS-G(I)»H(IG^l) 

IEXITs2 

IF(GS.GE.0.>G0T092 

CSOcGS 

ALPHAs-2,»DF/GS 

IF(ALPHA.GT.1.0) ALPHAsl.O 

DF=F 

TOTsO. 

30 CONTINUE 
lEXITsS 

IF<IFN.EG.MAXFM)G0TO92 

ICONsO 

lEXITsI 

DO 31 1=1,N 

Z=ALPHA»W{IS+I) 

IF(ABS(Z).GE.EPS<I))1C0N=1 

31 X(I)sX(I)^Z 

CALL FUNCUN.X.FY.G) 

IFNsIFN+1 

IF(IGAMMA.EG.O) GO TO 33 
DO 34 1x1,M 
2 xALPHA*W(IS*I) 

34 X(I)=X(I)-Z 
ALPHA=0.1«ALPHA 
PRINT 35.ALPHA 

35 FORMAT(1X,»ALPHA=*,F10.5) 

GO TO 30 

33 CONTINUE 
GTSsO. 

DO 32 -,N 

32 GYSxGyS4C(I)»NvIS+I) 
IF(FY.GE.F)G0T040 
IF(ABS(GYS/GS0).LE..9)G0T050 
IF(GYS.CT.0.)C0TO40 
TOTxTOT+ALPHA 

ZzlO. 

IFCGS.LT.GYS)ZxGYS/( CS-UYS) 
IF(Z.GT.10.)Zs10. 

ALPHAxALPHA*Z 
FxFY 
GSxCYS 
GOT030 

40 CONTINUE 

DO 41 1x1,N 

41 X(I)xX<I)-ALPHA*W(IS4.t) 
IF(IC0N.E0.0)C0T092 
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I 


2=3. • ( f-fy ) / alpha-*-gys4<;s 
ZZsSQRK Z»»2-GS»GYS) 

Zsi .-(GYS«.2Z-Z)/(2.»Z2+GYS-GS) 

ALPHA=ALPHA*Z 

GOT030 

50 CONTINUE 
ALPHA.TOT^ALPHA 
FaFY 

IF(ICON.E0.0)GOTO90 
DF=DF-F 
DGS=CYS-CSO 
DO 51 Is1,N 
W(IGC>I)=G(1) 

51 G(I)=-H<IG+I) 

IF( DGS4. ALPHA»GS0. GT. 0.) GOTO60 
C COMPLEMENTARY DFP FORMUU 
SIG=1./GSO 
1R=-IH 

CALL MC11A<H,N,G.SIG,H,IR,1,0.) 

DO 52 1=1, N 

52 G(I)aW(IGG+I)-W(IG+I) 
SIG=1./(ALPHA»DGS) 

IRs-IR 

CALL MCnA<H,N,G,SIO,W.IR,0,O.) 
C0T070 

60 CONTINUE 
C DFP FORMULA 

ZZ=ALPHA/(DGS-ALPHA*GS0) 

SIG=-ZZ 

CALL MCI1A(H,N,G.SIG.W,IR,1,IE-13) 
Z=DCS»ZZ-1. 

DO 61 1=1,N 

61 G(I)=W(IGG^I)^Z«W(IG*I) 

SIC=1./(ZZ«DGS«*2) 

CALL MC11A(H,N,G,SIG.H,IR.0,0.) 

70 CONTINUE 

DO 71 1=1,N 

71 G(I)=W(IGG*I) 

GOTO20 

92 CONTINUE 
DO 91 1=1,N 
91 G(I)»W(IG4l) 

90 CONTINUE 

IF(IPRINT.EQ.O)RETURN 
PRINT 1001,ITN,IFN,IEXIT 
PRINT 1002,F 
PRINT 1002,(X(I),1.1,N) 

PRINT 1002,(0(1),I.1,N) 

RETURN 

END 
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SUBROUTINE VFOIZ(N.X.PHI.GPHI) 
REAL X(1),GPHI(1) 

COHMON/VFO1C/F,M,K,IS,MK.NU 
COHMON/VF01D/G(50) 

COMMON/VF01E/C(150) 
COMMON/VF01F/GC(1250) 

COMMON/VFO1C/T(150) 
COMMON/IGAMMA/IGAMMA 
IF(MK.EQ.1)CALL VF01B (N.H.X) 
MK=1 

if(igamma.eq.i) return 

PHIsO. 

DO 10 1=1,N 

10 GPHKDsGd) 

DO 12 1=1,M 
CC=C(I)-T(I) 

IF(I.GT.K)CCsAMIN1(CC,OEO) 

Y=T(IS+I)»CC 

IF(Y.EQ.0E0)GOT012 

PHI=PHI-fY«CC 

DO 11 J=1,N 

11 GPHKJ)sGPHI(J)+Y»GC((1-1)«NU+J) 

12 CONTINUE 
PHIs.5E0»PHIt-F 
RETURN 

END 


76 




MC1U PACE 1 





SUBROUTINE MC11A{A,N,Z,SIC.W,IR.MK,EPS) 
DIMENSION A(1).Z(1),W<1) 

C UPDATE FACTORS GIVEN IN A BY SIG«Z«ZTRAKSPOSE 
IFCN.GT.DGOTOI 
IR=1 

A(1)=A(1)-*-SIG •Z(1)««2 
IF(A(1).GT.O.)RETURN 
A(1)=0. 

IRsO 
RETURN 
1 CONTINUE 
NPsN-*-1 

IF(SIG.GT.O.)GOTONO 
IF(SIG.EQ.O. .OR.IR.EQ.ORETURN 
TI=1./SIC 
IJ=1 

IF(MK.EQ.0)GOTO10 
DO 7 1=1 ,N 

IF(A(IJ).NE.0.)TIsTI+W(I)»»2/A(IJ) 

7 IJsIJ+NP-I 
G0TO20 

10 CONTINUE 

DO 11 Isl.N 

11 W<I)=Z(I) 

DO 15 1=1,N 
IPs 1+1 
V=W(I) 

IF(A(IJ).GT.0.)G0T012 

W<I)sO. 

IJ=IJ+NP-I 
GOTO15 

12 CONTINUE 
TI=TI+V**2/A<IJ) 

IF(I.EQ.N)G0T014 
DO 13 JsIP.N 
IJ=IJ+1 

13 W(J)=W(J)-V*A(IJ) 

IN IJ=IJ+1 

15 CONTINUE 

20 CONTINUE 
IFdR.LE.O )G0T021 
IF(TI.GT.0.)G0T022 
IF(HK-1>40,40,23 

21 TIsO. 

IR=-IR-1 

GOT023 

22 TI=EPS/SIG 
IF(EPS.EQ.0.)IR=IR-1 

23 CONTINUE 
MMfil 
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TIMsTI 
DO 30 1=1,N 
J=NP-I 
IJ=IJ-I 

IF(A(IJ).NE.O.)TIMxTI-W(J)••2/A(IJ) 

W(J)=T1 
30 TIsTIM 
G0T041 

40 CONTINUE 
MM=0 

TIM=1./SIG 

41 CONTINUE 
IJ=1 

DO 66 1=1 ,N 

IP=I+1 

V=Z(I) 

IF(A(IJ).GT.0.)G0T053 

IFdR.GT.O .OR.SIG.LT.O. .0R.V.EQ.0.)G0T052 
IR=1-IR 

A(IJ)=V*«2/TIM 
IF(I.EQ.N)RETURN 
DO 51 JsIP.N 
IJsIJ+1 

51 A(IJ)sZ(J)/V 
RETURN 

52 CONTINUE 
TI=TIM 
IJ=IJ+NP-I 
G0T066 

53 CONTINUE 
AL=V/A(IJ) 

IF(MM)54.54.55 

54 TI=TIM+V»AL 
G0T056 

55 TI=W(I) 

56 CONTINUE 
RsTI/TIM 
A(IJ)=A(IJ)*R 
IF(R.EQ.0.)G0T070 
IF(I.EQ.N)GOT070 
B*AL/TI 

IF(R.0T.4.)G0T062 
DO 61 JsIP.N 
IJ=IJ*1 

Z(J)=Z(J)-V«A(IJ) 

61 A(IJ)=A(IJ)+B»Z(J) 

GCT064 

62 GM.TIH/TI 

DO 63 J«IP,N 
IJ*IJ^1 
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y=A(ij) 

A(IJ)=B»Z(J)+Y«GM 
63 Z(.J;sZ(J)-V«Y 
6« CONTINUE 
riMsTI 
IJsIJ^I 
66 CONTINUE 
70 CONTINUE 

IF(IR.LT.O)IR=-IR 

RETURN 

C FACTORIZE A MATRIX GIVEN IN A 
ENTRY MCnB 
IR=N 

IF(N.GT.1)G0TO100 
IF(A(1).GT.0.)RETURN 
A(1)=0. 

IRsO 

RETURN 

100 CONTINUE 
NPsN+1 
Ilsl 

DO 10i» Is2,N 
AAsA(II) 

NIsIUNP-I 

IF(AA.GT.0.)G0T0101 

AdDsO. 

IRsIR-l 

IIsNI+1 

GOTOlOW 

101 CONTINUE 
IP=II+1 
IIsNl4.1 
JK=II 

DO 103 IJ=IP,NI 
V=A(IJ)/AA 
DO 102 IK=IJ,NI 
A(JK)sA<JK)-A(IK)»V 

102 JK=JK>*.1 

103 A(IJ)«V 

104 CONTINUE 
IF(A(II).GT.0.)RETURN 
AdDaO. 

IR»IR-1 

RETURN 

C MULTIPLY OUT THE FACTORS GIVEN IN A 
ENTRY MCI 1C 
IF(N.EO.I)RETURN 
NPoN+l 
IIxN*Np/2 
DO 202 HIPs2,N 
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JKsII 

NIsII.I 

II=II-NIP 

AAsA(II) 

IPrIUI 

IF<AA.OT.0.)GOTO2O3 
DO 204 IJsIP.NI 
204 A(lJ)aO. 

GOTO202 
203 COHTINUE 

DO 201 IJ=IP.NI 
V=A(IJ>»AA 
DO 200 IK=IJ.NI 
A<JK)=A(JK)*A(IK)»V 

200 JK=JK+1 

201 A(IJ)aV 

202 CONTINUE 
RETURN 

C MULTIPLY A VECTOR 2 BY THE FACTORS GIVEN IN A 
ENTRY MCI ID 
IF(M.GT.1)GOTO30O 
2(1)s2(1)*A(1) 

W(1)s2(1) 

RETURN 

300 CONTINUE 
NPsNfl 
Ilsl 
N1sN>1 

DO 303 1=1.N1 
Y=Z(I) 

IF(A(II).EQ.0.)COTO302 

IJsII 

IPxUI 

DO 301 J=IP,N 
IJ>IJ-t-1 

301 Y»Y+Z(J)»A(IJ) 

302 Z(I)sY«A(II) 

W(I)»Z(I) 

303 IIsII+NP-I 
Z(N)=Z(N)-A(II) 

W(N)sZ<N) 

DO 311 K=1,N1 
I=N-K 

II»II-HP*I 

IF(Z(I).EQ.O.)OOT0311 

IPxI^I 

IJ«ZI 

Y.Z(I) 

DO 310 J.IP.N 
IJ>IJ4.1 
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310 Z(J)=Z(J)^A(IJ)»2(I) 

311 CONTINUE 
RETURN 

C MULTIPLY A VECTOR 2 BY THE INVERSE OF THE FACTORS GIVEN IN A 
ENTRY MC1TE 
IF(IR.LT.N)RETURN 
W(1)=Z(1) 

IF<N.GT.1)GOTOilOO 

2(1)=2(1)/A{1) 

RETURN 
1(00 CONTINUE 

DO 1(02 Is2.N 
IJsI 
11 = 1-1 
V=Z(r) 

DO 1(01 Jsl.II 
V=V-A(IJ)»Z(J) 

1(01 IJ=IJ+N-J 
H(I)=V 
1(02 Z(I)sV 

Z(N)=2(N)/A(IJ) 

NPsN-fl 

DO 1(11 NIP=2,N 
IsNP-NIP 
II=IJ-NIP 
Vs2(I)/A(II) 

IP=I>1 

IJsII 

DO 1(13 JsIP.N 
IIsII+1 

mo V=V-A(II)»Z(J) 

411 Z(I)=V 
RETURN 

C COMPUTE THE INVERSE MATRIX FROM FACTORS GIVEN IN A 
ENTRY MCI IF 
IF<IR.LT.N)RETURN 
A(1)=1./A(1) 

IF(N.EQ.I)RETURN 
NP=N-f1 
NIsN-l 
11=2 

DO 511 1=2.N 
A(II)*-A(II) 

IJrIl4.1 

IF(I.EQ.N)C0TO502 

DO 501 JsI.NI 

IK«II 

JK>IJ 

V=A<IJ) 

DO 500 K»I,J 
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JK=JK+NP-K 

V=V^A(IK)»A(JK) 

500 IKsIK+1 
A(IJ)*-V 

501 IJsIJ+1 

502 CONTINUE 
A(IJ)=1./A(IJ) 

AA=A(IJ) 

IJsI 

IP=I*1 

NI=N-I 

DO 511 J=2,I 

V=A(IJ)*aA 

IK=IJ 

KsIJ-IPW 

I1=IJ-1 

nip=ni+i,i 

DO 510 JK=K,I1 
A(JK)zA(JK)-*-V«A(IK) 

510 IKsIK-fNIP-JK 
A(IJ)sV 

511 IJaIJ4.NP-J 
RETURN 
END 
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SUBROUTINE VEO<iA(N.A,IA,B.BL.BU.X.Q.LT.K.G) 

DIMENSION A( lA,1),B'1).BL(1),BU(1).X(1).LT( 1) ,G( 1) 

IS=N 

lASsN 

IVsN 

ICACsN-fN 
IDsICAC 
DO 9 1=1 .N 

9 G(I)s-B(I) 

DO 10 1=1,N 
X(I)sO. 

LT(I)=I 

G(ICAC+I)sA(I,I) 

IF(0..GE.BL(I).AND.O..LE.BU(I))G0TO10 
IF(O..LT.BL(I))X(I)=BL(I) 

IF(O..GT.BU(I))X(I)=BU(I) 

DO 12 J=1.1 

12 G(J)sG(J)+A(J,I)«X(I) 

IF(I-EQ.N)G0TO10 
Ilrl+I 

DO 11 J=II.N 
11 G(J)=G(J)+A(I.J)»X(I) 

10 CONTINUE 
KsO 

K1 = 1 

20 CONTINUE 
lOUTsO 
DELsO. 

DO 21 I=K1,N 
LI=LT(I) 

IF(X(LI).EQ.BL(LI).AND.G(I).GE.0.)GOTO21 
IF(X(LI).EQ.BU(LI).AND.G(I).LE.O.)G0T021 
IF(G(I).LT.0.)GOTO22 
ZsX(LI)-BL(LI) 

Jsl 

G0T023 

22 CONTINUE 
Z=B'J(LI)-X(LI) 

JsO 

23 CONTINUE 

IF(G(ICAC>I).LE.0.)G0T02« 

BETA=ABS(G(I))/G(ICAOI) 

IF(BETA.GE.Z)G0T024 

ZsBETA 

Ds.5*Z»ABS(G(I)) 

Js-1 
G0TO26 
2« CONTINUE 

D*Z«{ABS(G(I))-.5*Z»G(ICAOI)) 

26 CONTINUE 
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IF(D,LT.DEL)GOT021 

OEUD 

ALPHASZ 

lOUTsI 

IINsI 

IFCJ.LT.OIINsO 

LBsJ 

21 CONTINUE 

IF( lOUT. NE.0)G0T029 

27 CONTINUE 

0 = 0 . 

DO 28 1=1,N 
LIsLTd) 

28 Q=Q-*-X(LI)»(G(I).B(LI)) 
Q=.5«Q 

RETURN 

29 CONTINUE 
SIGsl. 

IF(G(I0UT).GT.0.)SIG=-1. 
LIOUT=LT(IOUT) 

LIINsLIOUT 
25 CONTINUE 

SASsGdCAC+IOUT) 
IF(K.EQ.O)GOT031 
DO 30 1=1,K 

30 GdS-t.I)sG(IEU.I)»AdOUT,I) 

31 CONTINUE 

DO 37 I=K1,N 
LIsLTd) 

IF(LI-LI0UT)32,37,33 

32 Z=A(LI,LI0UT) 

G0T034 

33 Z=A(LIOUT,LI) 

3^* CONTINUE 

IF(K.EQ.O)GOT036 
DO 35 J=1,K 

35 Z=Z-Ad,J)«GdS+J) 

36 GdS+DsZ 

37 CONTINUE 
GdSH.IOUT)=SAS 
IF(K.EQ.O)GOT0^2 
GdS+K)s-AdOUT,K) 
IF(K.EQ.1)C0T0i»2 
I=K 

DO U1 11=2,K 
1 = 1-1 

Z=-AdOUT,I) 

I1=I+1 

DO no JsIl.K 

no z=z-Gdr4.j)»A(j,i) 
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*<1 GdS+DsZ 
‘•a CONTINUE 

IF(SIG.EQ.1.)G0T051 

DO 50 1=1,H 

50 G(IS^I)=-G(IS^I) 

51 CONTINUE 
IF(K.EQ.0)G0T062 
DO 61 1=1,K 

IF(G(IS+I).EO.O.)GOT061 

LIsLTCI) 

Jsl 


ZsBL(LI)-X(LI) 
IF(G(IS+I) .LT.O.)GOT060 
JsO 


Z=BU(LI)-X(LI) 

60 CONTINUE 
ZsZ/GdS+I) 

IF( 2.GE.ALPHA)G0T061 

ALPHA=2 

LB=J 

IIM=I 


LIINsLI 

61 CONTINUE 

62 CONTINUE 

X(LIOUT)=X(LIOUT)-fSlG*ALPHA 

IF(K.EQ.O)GOT071 

DO 70 1=1,K 

LI=LT(I) 

70 X(LI)sX(LI)+ALPHA»G(lS+I) 

71 CONTINUE 


DO 72 I=K1,N 

72 G(I)=G(I)>ALPHA*G(IAS+I) 

IF(riN.EQ.0)G0TO90 
X(LIIN)=BL(LIIN) 
IF(LB.E0.0)X(LIIN)=BU(LIIN) 
IF<IIN.EQ.IOUT)GOTO20 
K2=K-1 


sg=g(id+iin; 

I1sIIN+1 
DO 80 IsIl.N 
80 G(IV+I)sA(I,IIN) 
IF(IIN.EQ.K)G0T086 
I2=IIN+2 
S0=1 ./SC 
DO 85 IsIIN,K2 
V=C(IV+I1) 
VDsV/G(ID^I 1) 
S1=S0-».V*VD 
RsSI/SO 

G(ID^I)=G(ID-^I1)*R 
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BETASVD/SI 
IF(R.GT.H.)G0T08«1 
DO 81 Jsl2.N 

81 G(lV4j)sG(IVW)-V»A(J,I1) 
IF(I1.GT.K2)G0r083 

DO 82 JsI1,K2 
JIsJ^I 

82 A{J.I)3A(J1.I1)4.BETA*G(mj1) 

83 CONTINUE 
A(K.I)sBETA 
DO 84 JsKI.N 

84 A(J.I)sA(J.I1)-i^BCTA«G(IV4.J) 

G0T0849 

841 CONTINUE 
IF(I1.GT.K2)G0T0843 
DO 842 JsI1,IC2 
JlsJ-fl 

842 A(J.I)sBETA*G(IV<*-J1)4>A(J1,I1)/R 

843 CONTINUE 
A(K.I)sBETA 
DO 844 JsKI N 

844 A(J.I)sBETA*G(IV 4-J)«A(J.I1)/R 
DO 845 JsI2.N 

845 G(IV+J)sG(IVW)-V»A<J,I1) 

849 CONTINUE 

LT(I)8LT(I1) 

SOsSI 

I1sI2 

85 I2sI 24>1 
SGsl./SI 
LT(K)sLIIN 
C(ID4 -K)sSG 
IF(IIN.EQ.1)G0T0851 
IIsIIN-1 

DO 852 I«1,II 
ZsA(IIN.I) 

DO 853 J=IIN,K2 
853 A(J,I)sA(J+l,I) 

852 A(K.I}>Z 
851 CONTINUE 

86 CONTINUE 

DO 87 laKI.N 

87 G(ICAC+I)sG(ICAOI)*SG»C(IV+I>«»2 
KlaK 

KsK2 

IINaO 

ALPHA.1E75 
SASsG(ICAOIOUT) 

IF(SAS.GT.O.)ALPHAsABS(G(lOUT))/SAS 
IF(G(lOUT).LT,0.)G0T0898 
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Ji1 

Z=X(LIOUT)-BL(LIOUT) 

G0T0899 

898 CONTINUE 
JsO 

Z=BU(LIOUT)-X(LIOUT) 

899 CONTINUE 

IF(Z.GE.ALPHA)G0T025 

ALPHASZ 

LBsJ 

IINsIOUT 

LIINsLIOUT 

G0T025 

90 CONTINUE 
K2sKU1 

IF(SIG.EQ.1.)G0T091 
DO 901 IsKI.N 
901 G(IAS+I)=-G(IAS>I) 

91 CONTINUE 
IF(I0UT.EQ.K1)G0T097 
LT(I0UT)sLT<K1) 

LT(K1)=LIOUT 

G(IAS+IOUT)sG(IAS«-K1) 

G( ICAC-^IOUT) sG( ICACU-K1) 

GdCACfKDsSAS 

G(IOUT)sG(K1) 

IF(K.EQ.0)GOTO97 

DO 92 1=1,K 
Z=A(K1,1) 

A(K1,I)=A(IOUT,I) 

92 AdOUT.DsZ 

93 CONTINUE 
IF(K2.EQ.I0UT)G0T095 
I1sIOUT-1 

DO 94 I=K2,I1 

94 AdOUT.DsACl.KD 

95 CONTINUE 

IF( I0UT.EQ.N)G0T097 

IlsIOUT^I 

DO 96 IsIl.N 

96 Ad,IOUT) = Ad.K1) 

97 CONTINUE 
G(K1)sO. 

KsKi 

IF(K.EQ.N)G0T027 
DO 98 I=K2,N 
ZsGdASfD/SAS 
Ad.KDsZ 

98 GdCAC+I)sGdCACi-I)-Z«C(IAS*I) 
K1sK2 
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GOTO20 

ENTRY VE. 3 

IF(K.EO.O>RETURN 

ID:N<».N 

G{N+1)sU/G(ID*1) 
IF(K.EQ.I)RETURN 
N1sK-1 
DO 111 
11 = 1+1 

A(I1.I)=-A(I1.I) 
IF(I.EQ.N1)GOT0102 
I1=1+2 

DO 101 J=I1.K 
Z=A(J.l) 

J1=J-1 

DO 100 LsIl.JI 

100 Z=Z+A(J,L)»A(L.I) 

101 A(J.I)=-Z 

102 CONTINUE 

AA=1./G(ID+I1) 

G(N+I1)sAA 

DO 111 Js1,I 

ZsA(I1,J)»AA 

G(N+J)sG(N+J)+Z»A(I1,J) 

IFd.EQ.DGOTOm 

JlsJ+1 

DO 110 LsJI.I 

110 A(L,J)=A(L,J)+A(I1.L)»Z 

111 A(I1,J)=Z 
RETURN 
END 
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PROGRAM MRRV(INPUT,OUTPUT.TAPE10) 

DIMENSION X(16),EPS(16) 

COMMON/NFEST/NFEST 
C0MM0H/VF01E/C(150) 

C0MM0N/VF01F/GC(25.50) 

C0MM0N/VF01G/T(150) 

C0MM0N/VF01I/G2P(325) 

C0HM0N/P/TF,NA,TTA(7).TA(7),NM,TTM(7> ,TM(7) 

COMMON/CNSTR/ME,Ml,MC 

COMMON/TIME/TMAX,TO, T1,IRS,MINS,AR,ITN,IFH 

COMMON/INOM/INOM 

COMMON/IPROS/IPROB 

DATA DPR/57.29577951/ 

IPR0B=2 

IPROBsl 

IRSsl 

IRSsO 

TMAX=1200. 

N=11 

NA=5 

NM:5 

M=1 

MCsM 

Ks1 

MEsK 

MIsM-K 

X(l)s3200./806.‘»8263 
DO 5 1=1,5 

TTA(I)srLOAT<I-l)»0.25 
X(I+1)=20./DPR 
TTM(I)=FL0AT(I-1)»0.25 
X(I+NA+l)s60./DPR 
5 CONTINUE 
DO 1 1=1,N 
EPS<I)=1.E-4»X(I) 

1 CONTINUE 
AKMINsl.E-4 
0FN=.5 
MAXFNs10000 
IPRlsl 
IPR2=1 
IWs2500 
H0DE=1 

C(M4.1)x0.136 

AKslE60 

M2=M+M 

NN=N»(N4.1)/2 

MINS=1 

ITNsO 

IFNsO 
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PRINT 20 
20 FORMAT(IHI) 

IF(IRS.EQ.O) GO TO 8 

READ(10-16) IRS.MINS.MODE,KFEST,ITH,IFH,(X(I),Is1,M),(C(I).1*1,H2 
1 ),(T(I).Is1,M2),(G2P(I).I.1.NM),AK.DFN 

16 FORMAT(9X,6(T10),135(/.9X.5020)) 

PRINT 17. IRS,MINS.NODE.NFEST,ITH.IFN,(X(I),1=1,H),(C(I),1*1,M2 

1 ).(TCI).1=1,M2).(02" T).I=1.MM).AK.DFN 

17 F0RMAT(9X.6(110).135(/.9X,5E2‘I.14>> 

PRINT 15 

15 FORMATC//IX*INPUT VALUES OF THETA AND SIGMA USED*) 

PRINT 22.NFEST 

22 F0RMAT(1X*INPUT VALUE OF HESSIAN USED*//1X* TOTAL FUNCTION • 

$"EVALUATIONS SO FAR*18) 

8 CONTINUE 

CALL SECOND(TO) 

CALL VF01A(N.M.K.X.EPS,AKHIN,DFM.MAXFM,IPR1,IPR2.IW.MODE) 

IF(T1-T0.LT,TMAX) GO TO 9 

IRS=1 

PRINT 30.TMAX 

30 F0RMAT(//1X».TIME LIMIT OF*G15.8*EXCEEDED.*//) 

PRINT 31,T1-T0 

31 F0RMAT(1X*C0MPUTATI0N TIME»C16.8/) 

MODEs-3 

9 REWIND 10 

WRITE(10.16) IRS.MINS.MODE.NFEST.ITN.IFN.(X(I) ,1=1 ,r.’) ,(C(1) ,1*1 .M2 
1 ).(T(I).I=1.M2),(G2P(I).I=1,MN),AK.DFN 

PRINT 17. IRS,MINS,MODE,NFEST,ITN.IFN.(X(I),1=1,N),(C(I),1=1.M2 

1 ),(T(I),1=1,M2),(G2P(I).1=1,NN).AK.DFN 

PRINT lO.NFEST 

10 F0RMAT(//5X*T0TAL FUNCTION EVALUATIONS*1X16) 

END 
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SUBROUTINE VFOIB(N.N.X) 

DIMENSION X(1) 

COMMON/NFEST/NFEST 

COMMON/VFO1C/F,HM,K,IS,MMK,MU 

COMMON/VF01D/G(50) 

C0MMOM/VF01E/C(15O) 

COMMON/VFO1F/GC(25.50) 

COMMON/IGAHHA/IGAHHA 
DATA NFEST/0/ 

CALL SG(X.F) 

IF(IGAHHA.EQ.1) RETURN 

NFEST=MFEST4l 

CALL SGX(N.X,M.NFEST) 

RETURN $ END 





I 


SG PAGE 1 




SUBROUTINE SG(X.F) 

REAL M.N.MASS 

dimension X(16).Y(9).2(9),DT(9) 
COMMON/VFO1E/S(150) 

COMMON/SS/SS(50) 
COMHON/INOM/INOM 
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C0MM0N/P/TF.NA.TTA(7).TA(7),HM,TTM(7),TM(7> 
COMMON/CNSTR/ME.MI.MC 


COMMON/IGX/IGX 

COMMON/TRAJ/A,M,N.KOUNT.TIMETB.MASS.NDE 
COMMON/COSI/CF.VCGOR.CSI,SSI 
COMMON/IPROB/IPROB 
COMMON/IGAMMA/IGAMMA 

DATA M,IGX,DPR/5.880964IE-2,0,57.29577951/ 

IGAMMAsO 

TF=X(1) 


DO 5 1=1.NA 
TA(1)=X(I>1) 


DO 6 I'l.NM 
6 TM(1)=X(I+NA*1) 
TIMETB=X(1+NA+NM+1) 

IF(IPR0B.EQ.2) GO TO 59 
TIMETS=10. 

59 CONTINUE 


N0Es6 

DELT=.004 $ NISs250 
TsO. 

T<1)=0.0 $ Y(2)a0.0 $ Y(3)=1.017432502 
Ym:o!T?(8)-J Y(5)=..0209«3951 $ Y(6)=0. 

KOUNTsO 
TIMETLsTIMETB 
DO 15 1=1,NIS 
IF(INOM.NE.I) go TO 25 
IF(M0D(I,25).NE.1) GO TO 25 
if(i.eq.i) print 60 

Ap- 6x :• 

39X,* (SEC)* ,3X,»(DEG)*,3X,*(DEG)*,3X,*(FT)* 3X •(FT/SECI* 

1 •(I>EG)..3X,A(DEG)-.3X,A(DE0)» 3X • 0EG)i) ’ ’ 

CALL DERIV(T,Y,DY) 

Z(1)=DPR*Y(1) $ Z(2)=DPR*Y(2) $ Z(3)=2.0926428E7*fYfi ^ 
z 4).259«7.772»X(A) « Z(5) = DP«.Y(5) » Z(i).DPm(6) 

Z(7)=806.48263*7(7) $ Z(8)=806.48263*7(8) 

PA=DPR*A $ PM=DPF*M 
TIME=806.48263*T*TF 
TEHPsVCCOR*CSI+W*CF 

CI=CF*TEMP/SQRT(TEMP**24.(VCG0R*SSI)*»2) 

PRINT Z0.T.TIM£,(Z(J).j=i,6),pa,PM.Z(7).Z(8).N 
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20 FORMAT(1X,F5.3.2X,F6.1,2X,F6.2.2X,F6.3.2X,F7.0,2X,F6.0,2X,F6.2,2X, 

1 F6.2,2X.F6.3,2X,F6.2.2X,F6.3.2X,F6.3.2X.F6.3> 

25 CONTINUE 

IF(KOUNT.EQ.2) GO TO 10 
DELT1=TIMETL/TF-T 
IF(DELTI.LE.DELT) GO TO 11 

10 CALL RUNGE(T.Y.DELT.NOE) 

IF(Y(5).GT.-1.) GO TO 50 
IGAMMAsI 

RETURN 
50 CONTINUE 
GO TO 15 

11 CONTINUE 

CALL RUNGE(T.Y.DELT1,NDE) 

KOUNT=KOUNT+1 
IF(INOH.NE.I) GO TO 63 
TIMEL=806.il8263»TIMETL 
IF(KOUNT.EQ.I) PRINT 61.TIMEL 
61 FORMATdH .•IGNITION OCCURRED AT •,F10.5,» SECONDS*) 

63 CONTINUE 

IF(DELTI.EQ.DELT) GO TO 12 

DELT2=DELT-DELT1 

CALL RUNGE<T,Y,DELT2,NDE) 

12 TIMETLsTIMETUO,63901697 
15 CONTINUE 

CALL DERIV<T.Y.DY> 

TEMPsVCGOR«CSI+W«CF 

Cl2CF»TEMP/SQRT(TEMP*»2+(VCG0R*SSI)»«2) 

IF(INOM.NE.I) GO TO 26 

Z(1) = DPR*Y(1) $ Z(2)sDPR*Y(2) $ Z(3)=2.0926'<28E7»(Y(3>-1.) 
Z(«)s25947.772»YCi«) $ Z(5)=DPR»Y(5) $ Z(6)=DPR»Y(6) 
Z{7)s806.48263*Y(7) $ Z(8)s806.«8263»Y(8) 

PAsDPR*A $ PMsDPR»M 
TIMEs806.48263»T»TF 

PRINT 20.T,TIME,(Z(J),J=1.6),PA.PM,Z(7),Z(8),N 
26 CONTINUE 
INOMsO 

IF(IPROB.EQ.2) GO TO 72 
F=-Y(2) 

SS(1)s(Y<3)-1.0047786464)/0.0047786464 
SS(2)=Y(7) 

SS(3)=Y(8) 

GO TO 74 
72 CONTINUE 
FsCI 

SS(1)=(Y(3)-1.029054170)/.029054170 
SS(2)=Y(4)/.90566542-1 . 

SS(3)=Y(5) 

.SS(4)=Y(7) 

SS(5)=Y(8) 
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CONTINUE 

if(igx.eq.i) go to iio 

DO 35 1=1,MC 
S(I)sSS(I) 

35 CONTINUE 


PRINT 28,F,(S<J) .Jsl.MC) 
28 F0RMAT(1X./,1X.6F20.15./) 
‘JO CONTINUE 


RETURN 

END 
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SUBROUTINE RUNGE (T.X.DELT.M) 

DIMENSION X(10),DX(10),DELX(10,3).XV(10) 

CALL DERIV (T,X,DX) 

T2 s T ♦ DELT/2.0 
DO 100 I s 1,N 
DEUd.l) = DX(I)»DELT 
100 XV(I) 8 X(I) ♦ DELX<I,1)/2.0 
CALL DERIV (T2,XV,DX) 

DO 200 Isl.N 
DELX(I,2) 8 DX(I) •DELT 
200 XV(I) 8 X(I) ♦ DELX<I.2)/2.0 
CALL DERIV (T2,XV,DX) 

DO 300 Isl,N 
DELX(I,3) = DX(I)»DELT 
300 XV(I) 8 X(I) ♦ DELX(I.3) 

T 8 T + DELT 

CALL DERIV (T.XV.DX) 

DO <400 Isl .N 

400 X(I) 3 X(I) ♦ (DELX(I.I) ■*. DX(I)»DELT ♦ 2.0»(DEU(I,2) ♦ DELX(I,3) 
1 ))/ 6.0 
RETURN 
END 
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SUBROUTINE DERIVCTT,Y,0Y) 

REAL L.H.N,NLIH,MASS 
DIMENSION Y(9).DY(9> 

COMMON/P/TTF.NA.TTA(7) .TA(7) ,1IM,ITM(7) ,TM(7> 

COMMON/TRAJ/A,M.N,KOUNT.TIHETB.MASS.NDE 
COMMON/COSI/CF,VCGOR.CS,SS 

DATA TW, W2 , ALIM, NLIM/. 11761928.3• 'I585739E-9..69813170,4.5/ 

TsY(1) $ FsY(2) $ R*Y(3) 8 V.Y(4> $ GsY(5) 8 SsY(6) 

SF=SIN(F) $ CFsCOS(F) 

SG>SIN(G1 8 CG:COS(G) 

SS:SIH(S) 8 CSsCOS(S) 

TFsSF/CF 8 TGsSG/CG 
CALL SLINKTT.A.TTA.TA.NA) 

CALL SLINKTT.M.TTM.TM.NM) 

SA:SIN(A) I CA:COS(A) 8 SMsSIN(H} 8 CHsCOS(M) 

VCG=V»CG 

VCGOR=VCG/R 

CALL AERG^R.V.A.SA.CA.D.L) 

GRsI./R»»2 

IF(KOUNT-Gr.O) GO TO 5 
MASSsI. 

THRsO. 

GO TO 15 

5 IF(KOUNT.GT.I) GO TO 10 

MASSs1,-.83533982«(TT«TTF-TIMETB) 

THRs.30555556 
GO TO 15 

10 MASSs.46620367 
THRsO. 

15 CONTINUE 
NsL/MASS 

TMDOMs(THR» CA-D)/MASS 

TPLOM=(THR*SAfL)/MASS 

TDOTsVCCOR*CS/CF 

FDOTsVCCOR«SS 

RDOTsV»SC 

VDOrsTMDOH-GR»SG+W2»R'»CF»(SG»CF-CG»SF*SS) 

GDOTs TPLOM* CM/ V-G R» CG/V+VCGOR+TW* CF» CS+W2* ( R» CF/V ) • ( CC» CF+SG^SF^SS 
1 ) 

SDOT=TPLOM«SM/VCG-VCGOR«TF»CS^TW»(TG»SS»CF-SF)-(W2/VCGOR)»SF*CF»CS 

DY(1 )=TDOT»TTF 8 DY(2)=FDOT*TTF 8 DY(3)sRDOT*TTF 

DY(«)=VDOTnTF 8 DY<5)=G[)OT*TTF 8 DY(6)sSDOT*TTF 

IF(HDE.EQ.6) GO TO 20 

Psi.-A/ALIM 

IFCP.GT.O.) PsO. 

DY(7)=-P*P 

PsM 

IFCP.GT.O-) P=0. 

DY(8)s-P«P 

DY(7) = DY(7)»TTF $ DY(8) = DY(8)nTF 
20 CONTINUE 

RETURN 8 END 
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SUBROUT INE AERO(R,V.AL,SAL,CAL.D.L) 

REAL L.M 

DIMENSION TM(5).TAA(5),TBB(5).TCC(5).TA1(5),TA2(5>.TB1(5),TB2(S),T 
AC1(5).TC2(5),TD1(5).TD2(5).TD3<5) 

COHHON/IPROB/IPROB 
DATA TM/.2.1.2.5..10.,20./ 

DATA TAA/0..0..0.,2.4E-12,1.38e-12/ 

DATA TBB/7.8E-8.7.7E-8.5.6E-8.-7.11E-7.-3.81E-7/ 

DATA TCC/.0114..0076..0028..0578..0297/ 

DATA TA1/-1.0e-4.-7.8lE-5.1.09E-5.1.41E-5.1.33E-5/ 

DATA TB1/-2.19E-3.-1.25E-3,-9.62E:-4,-9.0E-4,-8.19E-4/ 

DATA TA2/-1.0E-4.-7.8lE-5.4.86E-6,-6.6E.6.-6.25E-6/ 

DATA TB2/-2.19E-3.-1.25E-3,-1.13E-4,3.49E-4,3.58E-4/ 

DATA TCI/.035,.076,.0324,.0261,.0243/ 

DATA TC2/.035,.076,.0204,.0114,.0105/ 

DATA TD1/1.33E-4.-7.81E-5,2.94E-4,4.38E-4,4.38E-4/ 

DATA TD2/2.68E-2,2.90E-2,1.41E-2,6.50E-3,6.50E-3/ 

DATA TD3/-.030,-.030,-.035,-.035.-.035/ 

DATA DPR,S/57.29577951.125.84/ 

Hs20926428.«(R-1.) 

VEL=25947.772»V 
ALFA=DPR«AL 

CALL ATH62(H.TAU.SIGMA) 

MsVEL/< 65.770»SQRT(TAU)) 

CALL SLIN1(M,AA,TM,TAA,5) 

CALL SLIN1(M,BB,TM,TBB,5) 

CALL SLIN1(M,CC,TM,TCC,5) 

IF(ALFA.0T.16.)G0 TO 10 
CALL SLIN1(M,A,TM,TA1,5) 

CALL SLIN1(M,B.TM.TB1.5) 

CALL SLIN1(M,C,TM,TC1,5> 

GO TO 20 

10 CALL SLIN1(M,A,TM,TA2.5) 

CALL SLIN1(M,B,TM,TB2,5) 

CALL SLINUM,C,TM,TC2,5) 

20 CONTINUE 

CALL SLIN1(H,D1,TM,TD1,5) 

CALL SLIN1(M,D2,TM.TD2,5) 

CALL SLIN1(M,D3,TM,TD3,5) 

CASF=AA»H»H^BB«H4.CC 
CAPRs A» ALFA* ALFA+B* ALFA+C 
CArCASF<i>CAPR 

CNs D1•ALFA* ALFA+D2* ALFA+D3 
CL*CM*CAL-CA*SAL 
CDsCA*CAUCN*SAL 
IF(IPR0B.EQ.2) GO TO 30 
BBs15268.635*SIGMA*V**2 
GO TO 40 

30 BBs932.34367*SIGMA*V**2 

40 CONTINUE 
DsBB*CD 
LsBB*CL 
RETURN 
END 
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SUBROUTINE ATH62(H.T,S) 

REAL LH(8) 

DIMENSION HH(8).TH(8).CS(d) 

DATA HM/0.. 36089 • ,65617.. 10*1987., 15** 199., 17060*1. ,200131 • ,259186./ 
DATA TM/288.15,216.65,216.65.228.65,270.65,270.65.252.65,180.65/ 
DATA LM/-1.9812E-03.0..+3.0*I8E-04.8.53***»E-0*I,0.,-6.096E-0*I, 

1 -1.2192E-03,0./ 

DATA CS/3.401824655257E-11.1.683376997149E^00, 

19.817858914969E+80.1.506414967722E+29,4.39474988448IE-01, 
24.717851690435E-43,1.562793740651E-22,5.028946984109E+01/ 

DATA GOR/.010413309/ 

IFdl.GE.O.) GO TO 5 
T=288.15-1.9812E-3*H 
S=1.-2.926291 
RETURN 

5 IF(H.LE.750000.) GO TO 10 
T=180.65 
S:8.111816E-18 
RETURN 
10 CONTINUE 
DO 20 1:1.8 
20 IF(H.GE.HM(I))MsI 

TsTM(M)♦LM(M)•(H-HM(M)) 

IF(LM(M).EQ.0.)30.40 
30 SsCS(M)«EXP(-GOR»H/T) 

RETURN 

40 SsCS(M)«T»»(-1.-GOR/LM<M)) 

RETURN 

END 
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SUBROUTINE SLIHUX. Y. AX, AY. N> 

C SINGLE LINEAR INTERPOUTION SUBROUTINE 
DIMENSION AX(N), AY(N) 

IF(N.EQ.t) CO TO 30 
IFd.LE.AXd)) 00 TO 20 
IF(X.CT.AX(N))00 TO 30 
DO 10 Isl.N 
K « I 

Z ■ X - AX(1) 

IF(Z.LE.0.)OO TO 11 

10 CONTINUE 

11 J»K-1 

S s (AY(K) - AY(J))/(AX(K> - AX(a)) 

Y « AY(J) ♦ S«( X - AX(J)) 

GO TO 100 
20 YsAY(l) 

GO TO 100 
30 Y»AY(N) 

100 RETURN 
END 
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SUBROUTINE SGX(N.X.H.NFEST) 

DIMENSION X(16).XF(16),XB(16),CF(16),CB(16) 
COMM0N/VF01D/G(5O) 

COHMON/VFOlF/GC(25,50) 

COHMON/SS/SS(50) 

COMMON/IGX/IGX 
DATA EPSP/l.E-M/ 

IGXsl 

DO 15 Ia1.H 
XF(I)=X(I) 

15 XBCDsXd) 

DO 20 1=1.N 
DX=ABS(EPSP»X(I)) 

IF(ABS(X(I)).LE.EPSP) 0X=EPSP«2 
XF(I)sX(I)+DX 

CALL SG(XF.FF) 

DO 16 Jsl.M 

16 CF(J)=SS{J) 

XF(I)=X(I) 

XB(I)sX(I)-DX 
CALL SG(XB,FB) 

NFEST=MFEST4.2 
DO 17 Jsl.M 

17 CB(J)«SS(J) 

XBCDsXd) 

Gd)s0.5»(FF-FB)/DX 
DO 20 J=1,M 

GCd,J)s.5»(CF(J)-CB(J))/DX 
20 CONTINUE 
IGX=0 
RETURN 
END 




APPENDIX D 

GLOSSARY OF VARIABLES IN USER CODE 
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PR JGRAM MRRV 


Var l.iMi- 


Definition 

AK 

Real 

Measure of constraint convergence rate 

AKMIN 

Real 

Defines convergence for constraints 

C(I) 

Real 

Constraint residuals and constraint scale factors 

DFN 

Real 

Expected decrease In performance index 

DPR 

Real 

Degrees per radian 

EPS(I) 

Real 

Defines convergence in VA09A 

GC(I) 

Real 

Derivatives of the constraints 

G2P(I) 

Real 

Second derivative matrix in factored form 

I 

Integer 

Counter 

IFN 

Integer 

Number of function evaluations in VA09A 

INOM 

Integer 

Trajectory print flag. Print occurs if INOM • 1 

IPROB 

Integer 

Problem flag. ■ 1, reentry; • 2, plane change 

IPRl 

Integer 

Print flag in VFOIA 

IPR2 

Integer 

Print flag in VA09A 

IRS 

Integer 

Restart flag. * 0, first run; ~ 1, restart 

ITN 

Integer 

Number of iterations in VA09A 

IW 

Integer 

Number of elements in work space V(I) 

K 

Integer 

Number of equality constraints 

M 

Integer 

Number of constraints 

MAXFN 

Integer 

Maximum number of function evaluations In VA09A 

MC 

Integer 

- M 

ME 

Integer 

Number of equality constraints 

Ml 

Integer 

Number of Inequality constraints 


MINS 


Integer 


Number of iterations In V701A 


PROGRAM MRRV. Cont'd 


Variable 

Tjr£e 

Definition 

MODE 

Integer 

Flag which indicates how U, o, and G2P are input 

M2 

Integer 

- 2M 

N 

Integer 

Number of parametera 

MA 

Integer 

Number of parameters in angle of attack table 

NFEST 

Integer 

Number of function evaluations 

NM 

Integer 

Number of parameters in bank angle table 

NN 

Integer 

Number of elements In G2r 

T(I) 

Real 

Array where B and o are stored 

TA(I) 

Real 

Node ordinates for angle of attack table 

TP 

Real 

Pinal time 

TM(I) 

Real 

Node ordinates for roll angle table 

TMAX 

Real 

Internal time limit 

TTA(I) 

Real 

Node abscissas for angle of attack table 

TTM(l) 

Real 

Node abscissas for roll angle table 

TO 

Real 

Initial value of Internal time 

T1 

Real 

Current value of internal time 

X(I) 

Real 

Initial values of unknown parameters 

SUBROUTINE 

VPOIB 


Variable 

IZEl 

Definition 

C(I) 

Real 

Constraints and constraint scale factors 

F 

Real 

Performance index 

G(I) 

Real 

Derivatives of F with respect to X(I) 

GC(I,J) 

Real 

Derivatives of C(J) with respect to X(I) 

IGAMMA 

Integer 

Flag; ■ 0, Y > -60 deg ; ■ 1, y £ "60 deg 
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SUBROUTIHE VPOIB. Cont’d. 


Variable 

Type 

Definition 

IS 

Integer 

Used in optimization code 

K 

Integer 

Number of equality constraints 

M 

Integer 

Number of constraints 

MM 

Integer 

Used in optimization code 

MMK 

Integer 

Used in optimization code 

N 

Integer 

Number of parameters 

NFEST 

Integer 

Number of function evaluations 

NU 

Integer 

Used in optimization code 

X(I) 

Real 

Unknown parameters 

SUBROUTINE SG 

Variable 

Type 

Definition 

A 

Real 

Angle of attack (rad) 

CF 

Real 

Cos ^ 

Cl 

Real 

Cos 1 

CSI 

Real 

Cos Ip 

DELT 

Real 

Integration step size 

DELTl 

Real 

Partial integration step before engine is started 
or shut off 

DELT2 

Real 

Partial integration step after engine is started 
or shut off 

DPR 

Real 

Degrees per radian 

DY(I) 

Real 

Right-hand sides of the differential equations 
of motion 

F 

Real 

Performance index 
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SUBROUTINE SG. Cont'd. 


Variable 

Type 

Definition 

I 

Integer 

Counter 


IGAMMA 

Integer 

Flag; ■ 0, Y > -60 deg ; " 1, y ^ -60 deg 

IGX 

Integer 

Flag; - 1, prevents update of constraint residuals 
when derivatives are being computed 

INOM 

Integer 

Trajectory print flag; •> 1, causes trajectory to be 
printed 

IPROB 

Integer 

Problem flag; ■ 1, reentry; ■ 2, plane change 

J 

Integer 

Counter 

KOUNT 

Integer 

Engine flag; • 0, engine off; ~ 1, engine on; ■ 2, 
engine off again 

M 

Real 

Bank angle (rad) 

MASS 

Real 

Vehicle mass (nondlmensional) 

MC 

Integer 

Number of constraints 

ME 

Integer 

Number of equality constraints 

MI 

Integer 

Number of inequality constraints 

N 

Real 

Load factor 

NA 

Integer 

Number of parameters in angle of attack table 

NDE 

Integer 

Number of differential equations 

NIS 

Integer 

Number of integration steps 

NM 

Integer 

Number of parameters In bank angle table 

PA 

Real 

Angle of attack (deg) 

PM 

Real 

Bank angle (deg) 

S(I) 

Real 

Constraints (= C(I)) 

SS(I) 

Real 

Temporary values of constraints 

SSI 

Real 

Sin p 

T 

Real 

Time (normalized) 
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SUBROUTINE SG. 

Cont'd. 


Variable 

Type 

Definition 

TA(I) 

Real 

Node ordinates in angle of attack table (rad) 

TEMP 

Real 

Temporary computation 

TP 

Real 

Final time (nondlmenslonal) 

TIME 

Real 

Time (sec) 

TIMEL 

Real 

Ignition time or burnout time (sec) 

TIMETB 

Real 

Ignition tine (nondlmenslonal) 

TIMETL 

Real 

Ignition tine or burnout time (nondlmenslonal) 

TM(1) 

Real 

Node ordinates of bank angle table (rad) 

TTA(I) 

Real 

Node abscissas of angle of attack table (nondlmenslonal) 

TTH(I) 

Real 

Node abscissas of bank angle table (nondlmenslonal) 

VCGOR 

Real 

V cosy/r (nondlmenslonal) 

W 

Real 

Angular velocity of earth (nondlmenslonal) 

X(I) 

Real 

Unknown parameters 

Y(I) 

Real 

State variables (nondlmenslonal) 

Z(I) 

Real 

State variables (dimensional) 

SUBROUTINE RUNGE 


Variable 

Type 

Definition 

DELT 

Real 

Stepslze 

DELX(1,J) 

Real 

Increment in X(I) 

DX(I) 

Real 

Derivative of X(I) 

I 

Integer 

Counter 

M 

Integer 

Number of differential equations 
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SUBROUTINE RUNGE. Cont'd. 


Variable 

Type 

Definition 

T 

Real 

Independent variable 

T2 

Real 

Intermediate value of Independent variable 

X(l) 

Real 

Dependent variables 

XV(I) 

Real 

Intermediate value of X(I) 

SUBROUTINE 

DERIV 


Variable 

Type 

Definition 

A 

Real 

Angle of attack (rad) 

ALIM 

Real 

Angle of attack limit (rad) 

CA 

Real 

Cos o 

CF 

Real 

Cos 

CG 

Real 

Cos Y 

CM 

Real 

Cos p 

CS 

Real 

Cos 

D 

Real 

Drag (nondlmenslonal) 

DY(I) 

Real 

Right-hand sides of the equations of motion 

F 

Real 

Latitude, ^ (rad) 

FDOT 

Real 

e 

(nondlmenslonal) 

GR 

Real 

Acceleration of gravity (nondlmenslonal) 

KOUNT 

Real 

Engine flag; - 0, off; ■ 1, on; ■ 2, off again 

L 

Real 

Lift (nondlmenslonal) 

M 

Real 

Bank angle (rad) 



3 

S 
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SUBROUTINE DERIV, Cont'd. 


Vn r 1 /lit 1 1 ' 


UefInitInn 

MASS 

Real 

Mass (nondinenslonal) 

N 

Real 

Load factor 

NA 

Integer 

Number of nodes In angle of attack table 

NDE 

Integer 

Nuntber of differential equations 

NLIM 

Real 

Load factor limit 

NM 

Integer 

Number of nodes In bank angle table 

P 

Real 

Temporary computation 

R 

Real 

Distance from center of earth (nondlinenslonal) 

RDOT 

Real 

r (nondlmanslonal) 

S 

Real 

Heading angle, )|/ (rad) 

SA 

Real 

Sin a 

SDOT 

Real 

• 

<1/ (nondimenslonal) 

SF 

Real 

Sin 4> 

SG 

Real 

Sin Y 

SM 

Real 

Sin u 

ss 

Real 

Sin i|; 

T 

Real 

Longitude, 9 (rad) 

TA(I) 

Real 

Node ordinates in angle of attack table (rad) 

TOOT 

Real 

9 (nondlmenslonal) 

TP 

Real 

Tan ^ 

TG 

Real 

Tan Y 

THR 

Real 

Thrust (nondlmenslonal) 

TIMETB 

Real 

Ignition time (nondlmenslonal) 

TM(I) 

Real 

Node ordinates In bank angle table (rad) 
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SUBROUTINE DERIV, Cont'd. 


Variable 


Definition 

■ri’i.oM 

Kf.i 1 

(T Bln ti + L)/m 

TMDOM 

Real 

(T coa a - D)/m 

TT 

Real 

Normalized time, t 

TTA(I) 

Real 

Node abaclssas in angle of attack cable 

TTF 

Real 

Final time (nondimenaional) 

TTM(I) 

Real 

Node abscissas In bank angle table 

TW 

Real 

2ii) (nondlmensional) 

V 

Real 

Velocity (nondlmenslonal) 

VCG 

Real 

V cos Y 

VCGOR 

Real 

V coa y/r 

VDOT 

Real 

V (nondlmenslonal) 

W2 

Real 

( 0 ^ (nondlmenslonal) 

Y(I) 

Real 

Nondlmenslonal states 

SUBROUTINE 

_^R0 


Variable 

lype 

Definition 

AL 

Real 

Angle of attack (rad) 

ALFA 

Real 

Angle of attack (deg) 

A1 

Real 

Coefficient In C. 



SF 

A2 

Real 

Coefficient In C. 

A3 

Real 

rK 

Coefficient In C„ 

N 

B1 

Real 

Coefficient In C 



A. _ 

SF 

P2 

Real 

Coefficient in C. 

B3 

Real 

rK 

Coefficient In C„ 

N 

CA 

Real 

Axial force coefficient 
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SUBROUTINE AERO. 

Cont’d. 


Vdt I <it> 1 •• 

■'yp*- 

Ocf fnt cjfm 

CAL 

Real 

Cob a 

CAPR 

Real 

Axial force coefficient, preseure 

CASF 

Real 

Axial force coefficient, skin friction 

CD 

Real 

Drag coefficient 

a 

Real 

Lift coefficient 

CN 

Real 

Normal force coefficient 

Cl 

Real 

Coefficient in C, 

sr 

C2 

Real 

oefficleiit in CAp^ 

C3 

Real 

Coefficient In Cjj 

D 

Real 

Drag (nond linen 8Iona 1) 

DPR 

Real 

Degrees per radian 

H 

Real 

Altitude (ft) 

IPROB 

Integer 

Problem flag; - 1, reentry; • 2, plane change 

L 

Real 

Lift (nondimenBional) 

M 

Real 

Mach nund^er 

R 

Real 

Distance from center of earth (nondlmenslonal) 

SAL 

Real 

Sin a 

SIGMA 

Real 

Density ratio 

TAU 

Real 

Absolute temperature (deg K) 

TAl(I) 

Real 

Table for A1 

TA3(I) 

Real 

Table for A3 

TBl(I) 

Real 

Table for Bl 

TB3<I) 

Real 

Table for R3 

TCI(I) 

Real 

Table for Cl 
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SUBROUT INE^ A_ERO, Con t ' d . 


Variable 

T.yp.g 

TC3 

Real 

TM(I) 

Real 

T1A2(I) 

Real 

TxB2(I) 

Real 

T1C2(I) 

Real 

T2A2(I) 

Real 

T2B2(I) 

Real 

T2C2(I) 

Real 

V2L 

Real 

V 

Real 


Definition 

Table for C3 
Table for Mach nuaber 
Table for A2, a ^ 16 deg 
Table for B2, a ^ 16 deg 
Table for C2, a ^ 16 deg 
Table for A2, a > 16 deg 
Table for B2, a > 16 deg 
Table for C2, a > 16 deg 
Velocity (ft/sec) 

Velocity (nondimenslonal) 


SUBROUTINE ATM62 



Variable 

Type 

Definition 

CS(I) 

Real 

Coefficients lii density formula (nondimenslonal) ‘ 

CX)R 

Real 


H 

Real 

Altitude (ft) 

HM(I) 

Real 

Altitude at beginning of layer m (ft) 

I 

Integer 

Counter : 

Temperature gradient in layer m (deg K/ft) ^ 

LM 

Real 

M 

Real 

Denotes layer 

S 

Real 

Density ratio 

T 

Real 

Absolute temperature (deg R) 

TM 

Real 

Temperature at beginning of each layer (deg K) | 
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SIIBROITTINE 

SLINl 


Variable 

Type 

Definition 

AX(I) 

Real. 

Node abaciasaa 

AY (I) 

Real 

Node ordinates 

1 

Integer 

Counter 

J 

Integer 

Counter 

K 

Integer 

Counter 

N 

Integer 

Number of points in table 

S 

Real 

Slope 

X 

Real 

Abscissa 

Y 

Real 

Ordinate 

Z 

Real 

Temporary computation 
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